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A.  SCIENTIFIC  OBJECTIVES:  The  objective  of  this  research  program  is  to  devise 
innovative  joint  time-frequency  (JTF)  processing  concepts  for  radar  image  enhancement 
and  physics-based  feature  extraction.  In  particular,  we  shall  investigate  how  JTF 
techniques  can  be  utilized  to  enhance  synthetic  aperture  radar  (SAR)  and  inverse 
synthetic  aperture  radar  (ISAR)  imageries  by  removing  artifacts  due  to  uncompensated 
target  motion,  complex  target  scattering  physics,  articulating  target  components  and 
clutter  and  propagation  effects.  Furthermore,  we  set  out  to  re-interpret  the  extracted 
artifacts  in  a  more  meaningful  feature  space  so  that  they  can  be  utilized  to  enhance  the 
performance  of  target  identification  algorithms.  This  research  is  leveraged  against  our 
previous  JTF  work  under  the  Joint  Services  Electronics  Program,  as  well  as  our  state-of- 
the-art  radar  signature  simulation  capabilities.  At  the  end  of  this  program,  we  plan  to 
develop  a  set  of  JTF-based  radar  image  processing  tools  in  the  standard  MATLAB 
environment,  so  that  this  work  can  be  easily  disseminated  to  the  radar  community. 

B.  SUMMARY  OF  RESULTS  AND  SIGNIFICANT  ACCOMPLISHMENTS: 

During  the  first  year  of  this  research  program,  three  areas  of  investigation  have  been 
initiated.  First,  we  have  devised  a  joint  time-frequency  technique  to  remove  the 
interference  due  to  fast  rotating  parts  from  the  ISAR  imager  of  a  target.  We  utilize  a 
chirp-based  JTF  algorithm  to  extract  Doppler  lines  in  the  ISAR  image  and  estimate  the 


spin  rate  of  the  rotating  part.  Second,  we  have  processed  the  measurement  data  from  the 
TIRA  radar.  We  apply  the  previously  developed  adaptive  JTF  algorithm  to  form  high- 
resolution  ISAR  images  of  several  air  targets.  Third,  we  have  begun  a  study  on  using  the 
adaptive  wavelet  packet  transform  for  clutter  suppression  in  SAR  imagery.  We  show  that 
an  improved  target  signal-to-clutter  ratio  can  be  achieved  in  the  wavelet  domain.  The 
detailed  descriptions  of  these  three  projects  are  described  below. 

Removing  Image  Artifacts  due  to  Articulating  Target  Components.  It  is  well  known 
that  when  rotating  components  exist  on  a  target  such  as  gimbaled  antennas  or  propeller 
blades,  image  artifacts  are  introduced  in  the  Doppler  dimension  of  the  ISAR  image  [1,2]. 
These  smeared  features  oftentimes  overshadow  the  target  geometrical  features  and  hinder 
the  proper  interpretation  of  the  ISAR  image.  We  have  developed  a  technique  to  remove 
such  Doppler  smear  and  produce  a  clear  ISAR  image  of  the  target  based  on  adaptive  joint 
time- frequency  processing  [3,  4],  The  technique  entails  adaptively  searching  for  the 
linear  chirp  bases  that  best  represent  the  time-frequency  behavior  of  the  signal.  This  is 
accomplished  by  projecting  the  signal  onto  all  possible  chirp  bases  and  finding  the  one 
with  the  maximum  projection  value,  After  the  optimal  basis  is  found,  the  signal 
component  associated  with  this  basis  is  subtracted  from  the  original  signal.  By  iterating 
this  search  procedure,  the  signal  can  be  fully  parameterized  with  a  set  of  chirp  basis 
functions.  Since  the  Doppler  frequency  due  to  the  rotating  component  is  both  larger  and 
more  rapidly  varying  (in  dwell  time)  than  that  from  the  target  body,  the  signal 
components  due  to  the  fast  rotating  part  are  associated  with  those  chirp  bases  having 
large  displacement  and  slope  parameters.  On  the  other  hand,  the  signal  components  due 
to  the  target  body  motion  are  represented  by  those  chirp  bases  with  relatively  small 
displacement  and  slope  parameters.  By  sorting  these  chirp  bases  according  to  their  slopes 
and  displacements,  the  scattering  due  to  the  fast  rotating  part  can  be  separated  from  that 
due  to  the  target  body.  Consequently  a  cleaned  ISAR  image  of  the  target  body  can  be 
reconstructed  by  using  only  those  bases  associated  with  the  body  motion.  Furthermore, 
robust  Doppler  information  extraction  can  be  achieved  by  applying  the  period  detection 
algorithm  to  the  component  associated  with  rotating  parts  only.  We  have  applied  this 
algorithm  to  both  simulated  data  from  the  radar  scattering  prediction  code  Xpatch  [11] 


autocorrelation  function  versus  dwell  time  of  the  extracted  blade  contribution.  From  this 
plot,  we  can  determine  the  periodicity  of  the  signal  to  be  0.056  s.  The  reciprocal  of  the 
period  is  17.7  rps  and  it  corresponds  to  the  product  of  the  blade  rotation  rate  and  the 
number  of  blades  for  this  target.  We  believe  this  algorithm  can  also  be  extended  to  deal 
with  jet  engine  modulation  (JEM),  although  the  jet  problem  is  more  challenging  due  to 
the  high  rotation  rate  of  the  blades  and  the  electromagnetic  propagation  effect  through  the 
inlet  duct  structure.  This  topic  is  currently  under  investigation. 

Application  of  JTF  Motion  Compensation  Algorithm  to  NATO  Data.  We  have 
devoted  some  efforts  to  process  the  TIRA  radar  data  taken  in  Germany  in  November 
1997.  The  objectives  of  this  effort  are  to  test  our  previously  developed  adaptive  joint 
time-frequency  (AJTF)  algorithm  [10]  for  ISAR  motion  compensation  and  to  identify 
needed  areas  of  research  in  ISAR-based  target  recognition.  Our  algorithm  uses  a  search 
and  projection  technique  in  the  joint  (dwell  time)-(Doppler  frequency)  plane  to  select  and 
track  the  prominent  point  scatterers.  The  higher-order  translation  and  rotation  motions 
are  then  extracted  and  compensated  for  in  the  data  to  form  a  focused  image  of  the  target. 
The  motion  compensated  images  are  compared  against  the  reference  images  generated  by 
using  the  motion  information  available  in  the  instrumented  ARDS  data.  Furthermore, 
comparison  is  also  made  with  the  simulated  ISAR  images  of  the  air  target  from  the  radar 
signature  prediction  code  Xpatch.  A  set  of  representative  images  is  shown  below.  Figs. 
2(a),  (b)  and  (c)  are  the  images  for  a  target  at  azimuth=56°  (0°  being  nose-on)  generated 
using,  respectively,  AJTF  motion  compensation,  ARDS  sensor  information  and  Xpatch 
simulation.  The  dynamic  range  of  the  displayed  images  is  55  dB.  The  look  angle 
information  was  deciphered  from  the  ARDS  data.  By  comparing  Figs.  2(a)  and  2(b),  we 
observe  that  the  motion  compensated  image  and  the  ARDS -derived  reference  image 
appear  to  be  in  good  agreement.  While  this  is  true  for  most  angles,  we  did  find  that  the 
performance  of  our  motion  compensation  algorithm  depends  on  the  availability  of  a  clear 
point  scatterer  in  the  selected  range  cell.  A  general  criterion  should  be  further  developed 
to  automatically  select  the  range  cell  that  contains  a  strong,  well-isolated  point  scatterer. 
This  topic  is  being  further  pursued.  When  the  Xpatch  simulation  is  compared  against  the 
ARDS-derived  reference  image,  the  agreement  in  the  prominent  target  features  is  fair. 
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Fig.  2.  A  TIRA  target  at  azimuth=56°. 


(c)  Synthetic  image  from  Xpatch 
simulation. 


The  Xpatch  image  is  more  focused  and  does  not  exhibit  the  diffused  characteristics  of  the 
measurement  data.  A  more  detailed  CAD  model  should  improve  the  quality  of  the 
predicted  image.  Furthermore,  the  measured  images  contain  strong  jet  engine  modulation 
(JEM)  lines  in  the  frontal  look  region  of  the  aircraft.  Since  these  features  are  sensitive  to 
the  engine  spin  rate  relative  to  the  radar  pulse  repetition  frequency,  an  accurate  prediction 
is  difficult  by  Xpatch  simulation.  We  plan  to  use  this  set  of  processed  TIRA  data  to  test 
our  algorithm,  to  be  developed  under  topic  1,  for  removing  JEM  artifacts  from  the 
measured  data  in  order  to  unmask  the  geometrical  features  of  the  target. 


Clutter  Reduction  for  Radar  Images  Using  Adaptive  Wavelet  Packet  Transform. 

Synthetic  Aperture  Radar  (SAR)  images  of  ground  targets  generally  consist  of  target 
features  and  clutters  from  background  scattering.  In  automatic  target  recognition  (ATR) 
applications,  it  is  desirable  to  remove  the  clutter  from  the  target  image  before  ATR 
processing.  The  standard  way  to  suppress  clutter  is  to  apply  an  appropriate  threshold  level 
to  the  whole  SAR  image.  However,  this  approach  assumes  that  the  target  signal-to- 
clutter  ratio  (SCR)  is  large  enough.  Otherwise  this  direct  threshold  approach  results  in 
either  target  feature  loss  or  remnant  clutter  residue.  In  this  work,  we  set  out  to  develop  a 
decluttering  algorithm  to  automatically  extract  the  target  image  from  a  SAR  image  by 
maximizing  the  SCR  using  the  adaptive  wavelet  packet  transform  (AWPT)  [5].  The 
wavelet  packet  basis  is  a  generalization  of  the  conventional  wavelet  basis  [6]  and  has 
been  applied  for  image  compression  [7]  and  moment  matrix  sparsification  [8].  Our 
approach  is  to  transform  the  SAR  image  to  a  new  domain  using  the  wavelet  packet  basis. 
Since  a  typical  target  image  usually  consists  of  point  scatterers  and  more  diffused  region 
features,  the  multi-scaled  wavelet  basis  is  well  suited  to  focus  the  target  image.  Clutter 
image,  on  the  other  hand,  is  statistically  uncorrelated  from  pixel  to  pixel,  and  the 
transformed  clutter  image  under  the  same  set  of  bases  remains  unfocused.  Therefore,  we 
expect  that  the  SCR  can  be  increased  by  transforming  the  original  image  using  an 
appropriately  chosen  set  of  wavelet  packet  basis.  The  cost  function  of  our  AWPT 
algorithm  is  chosen  to  describe  how  well  the  target  signal  is  focused  in  the  transform 
domain.  An  efficient  basis  search  algorithm  is  implemented  to  find  the  best  wavelet 
packet  basis.  Our  algorithm  is  tested  using  the  MSTAR  SAR  data  set  [9].  Fig.  3(a) 
shows  an  MSTAR  image  in  which  the  target  is  a  ground  vehicle  and  the  clutter  is  due  to 
vegetation.  There  are  several  strong  point  scatters  in  the  front  of  the  vehicle,  but  the 
scattering  from  the  back  part  is  relatively  weak.  Fig.  3(b)  shows  the  result  of  applying 
the  direct  thresholding  method  to  the  image.  Fig.  3(c)  shows  the  decluttered  image  by 
applying  the  AWPT  algorithm.  We  choose  Daubechies  filter  with  order  of  6  as  the 
wavelet  filter.  By  visually  comparing  Figs.  3(b)  and  3(c),  we  note  that  some  crucial 
features  of  the  target  are  kept  in  the  AWPT-processed  image.  In  both  processing 
methods,  there  is  some  target  information  loss.  Fig.  3(d)  shows  the  signal-to-clutter  ratio 


Fig.  3(a).  SAR  image  of  a  ground 
vehicle  with  clutter 


Fig.  3(b).  Clutter  rejection  using 
the  direct  thresholding  method. 


Fig.  3(c).  Clutter  rejection  using 
the  AWPT  approach. 


Fig.  3(d).  SCR  vs.  target  image  loss  for 
the  two  processing  methods. 


versus  average  target  image  loss  for  the  two  processing  methods.  It  is  observed  that  for  a 
fixed  target  image  loss  the  AWPT  method  always  achieves  a  higher  SCR  value  than  the 
direct  thresholding  method.  Similar  results  are  obtained  when  the  algorithm  is  applied  to 
other  MSTAR  targets.  We  plan  to  examine  SAR  data  from  other  available  sources,  as 
well  as  to  examine  the  applicability  of  this  technique  to  ISAR  data. 


C.  FOLLOW-UP  STATEMENT: 


During  the  coming  year,  our  research  efforts  will  be  devoted  to  three  areas.  First, 
we  will  continue  our  work  on  extracting  Doppler  effects  in  ISAR  images.  We  plan  to 
develop  algorithms  for  the  more  challenging  problem  of  jet  engine  modulation.  The 
algorithms  will  be  tested  using  the  processed  TIRA  images  that  we  have  generated  during 
the  first  year.  Second,  we  will  devise  improved  algorithms  for  ISAR  motion 
compensation  to  overcome  computational  time  and  robustness  issues.  We  also  plan  to 
develop  a  variable  dwell  time  algorithm  to  maximize  cross  range  resolution  and  combat 
target  motion  jitters  encountered  in  measured  ISAR  data.  It  is  anticipated  that  a  new 
measurement  data  set  from  the  MERIC  radar  will  be  made  available  to  us  for  processing 
and  testing.  Third,  we  will  continue  our  initial  study  to  address  the  clutter  suppression 
issue  using  joint  time-frequency  concepts.  We  will  begin  to  examine  ISAR  data 
containing  clutter  by  requesting  and  processing  ship  data  from  the  Navy  small  craft  ATR 
program. 
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ISAR  Motion  Compensation  Via  Adaptive  Joint 
Time-Frequency  Technique 


A  novel  approach  for  inverse  synthetic  aperture  radar  (ISAR) 
imaging  is  presented  for  both  target  translational  motion  and 
rotational  motion  nonuniformity  compensation.  The  basic  idea 
is  to  perform  Doppler  tracking  to  individual  scatterers  via  an 
adaptive  joint  time-frequency  (AJTF)  projection  technique.  After 
maximizing  the  projection  of  the  phase  function  to  a  set  of  basis 
functions  in  time-frequency  plane,  the  Doppler  frequency  drift  of 
the  strongest  scatterer  in  the  range  bin  is  automatically  tracked 
out  and  the  multiple  prominent  point  processing  (PPP)  scheme  is 
implemented  to  eliminate  both  the  translational  motion  error  and 
rotational  motion  nonuniformity.  Further  the  azimuth  spacing  can 
be  estimated,  which  permits  polar  reformatting  of  the  original 
collected  data. 


I.  INTRODUCTION 

In  real-world  ISAR  (inverse  synthetic  aperture 
radar)  imaging  scenarios,  the  target  being  imaged  is 
often  engaged  in  complicated  maneuvers  that  combine 
translational  and  rotational  motions.  Unless  a  good 
motion  compensation  algorithm  is  implemented, 
serious  blurring  can  result  in  the  ISAR  image  formed 
by  Fourier  processing.  A  motion  compensation 
algorithm  usually  consists  of  two  parts,  range 
alignment  and  cross  range  (Doppler)  tracking.  For 
various  existing  motion  compensation  algorithms, 
the  manner  in  which  range  alignment  is  performed 
is  fairly  standard.  It  is  accomplished  by  tracking  the 
time  history  of  a  reference  point  (such  as  the  peak  or 
the  centroid)  in  the  range  compressed  data  and  fitting 
it  to  a  polynomial  [1].  There  are,  however,  many 
different  schemes  to  perform  Doppler  tracking,  such 
as  the  subaperture  approach  [2-4],  the  cross-range 
centroid  tracking  approach  [5]  and  the  phase  gradient 
autofocus  technique  (PGA)  [6].  All  of  these  methods 
consider  the  Doppler  frequency  shift  for  the  target 
as  a  whole,  and  apply  the  same  correction  vector 
to  all  the  scatterers  in  the  image.  However,  when 
the  dwell  time  is  long  or  when  the  target  exhibits 
fast  maneuvers,  the  phase  error  will  not  be  the 
same  for  all  the  scatterers  in  the  image.  Therefore,  a 
technique  which  can  automatically  track  the  Doppler 
frequency  shifts  of  individual  scatterers  is  needed 
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i  acquire  the  complete  information  about  the  target 
lotion.  We  propose  a  technique  which  utilitzes  joint 
me-frequency  concepts  to  achieve  this  goal. 

Joint  time-frequency  techniques  have  previously 
een  used  to  address  the  motion  compensation  issue 
7],  By  applying  the  time-frequency  distribution 
dries  (TFDS)  [8]  in  place  of  the  Fourier  transform 
ngine,  the  ISAR  image  can  be  effectively  examined 
t  each  time  instant,  thus  eliminating  range  drift 
nd  cross-range  smearing.  Unfortunately,  the  TFDS 
dchnique  is  based  on  the  Wigner-Ville  distribution 
nd  therefore  does  not  preserve  the  phase  information 
>f  the  original  image.  In  addition,  the  cross-range 
esolution  achievable  in  this  manner  is  still  less  than 
he  full  angular  aperture  of  the  original  data.  Recently, 
he  adaptive  joint  time-frequency  (AJTF)  concept 
vas  introduced  by  [9  and  10].  One  application  of 
his  concept  in  radar  image  processing  is  to  separate 
he  scattering  mechanisms  and  resonance  mechnisms 
if  ISAR  image  [11].  The  basic  idea  of  AJTF  is 
idaptively  looking  for  the  basis  functions  which 
lest  represent  the  joint  time-frequency  behavior 
>f  the  signal,  thus  parametrizing  the  signal  using 
hese  bases.  Here  the  AJTF  technique  is  used  as  a 
powerful  tool  to  perform  Doppler  frequency  tracking 
[12].  By  using  a  search  and  projection  procedure  in 
:he  time-frequency  plane,  the  reference  points  can 
automatically  be  selected  and  the  desired  motion 
parameters  can  be  figured  out,  thus  all  the  motion 
error  can  be  well  eliminated  by  multiplying  the  phase 
correction  vector,  interpolating  to  the  uniform  azimuth 
scale,  and  polar  reformatting  the  original  collected 
data.  The  original  phase  information  of  the  image  is 
kept  and  full  aperture  resolution  is  achieved. 

This  paper  is  organized  as  follows.  In  Section  II, 
we  introduce  the  mathmatical  model  of  motion  errors, 
[f  the  target  only  contains  translational  motion,  the 
phase  error  can  be  well  compensated  by  tracking 
one  reference  point  using  the  AJTF  technique  in 
she  target.  Further,  if  the  target  also  includes  the 
nonuniform  rotational  motion,  using  AJTF  technique 
so  track  two  or  three  different  reference  points,  a  more 
complex  motion  compensation  model  like  multiple 
orominent  point  processing  (PPP)  algorithm  [13,  14] 
can  be  utilized  to  compensate  both  translational  and 
rotational  motion  error.  But  unlike  what  original  PPP 
requires,  here  no  multiple  dominant,  well-isolated 
coint  scatterers  in  their  respective  range  cells  are 
needed.  Since  the  azimuth  scale  factor  can  also  be 
obtained  in  this  way,  the  polar  reformatting  is  possible 
or  the  imaging  of  large  targets  [2]. 

I.  ANALYTICAL  MODEL 

We  assume  that  the  target  rotation  outside  of 
he  plane  between  the  radar  and  target  is  small,  and 
tandard  range  alignment  has  been  applied  to  the 
lata  so  that  the  coarse  translational  motion  has  been 


removed.  All  the  scatterers  stay  in  the  right  range 
cells.  Then  the  phase  variation  as  a  function  of  dwell 
time  t  (time  index  of  each  pulse)  in  a  particular  range 
cell  x  can  be  written  as 

A4 

f{x,t)  =  YlAk 
k  =  1 

x  exp  \-j^^-(R(t)  +  xkcos6(0  +  yt  sin  0(f)) 

(1) 


where  Nk  is  the  number  of  point  scatterers  and  Ak  is 
the  magnitude  of  the  ifcth  point  scatterer.  R(t)  is  the 
residual  uncompensated  translational  displacement 
and  0(f)  is  the  rotational  displacement.  Note  that  either 
R(t)  or  0(f)  is  no  longer  the  simple  constant  or  linear 
relationship  which  is  assumed  for  ideal  no-motion 
error  ISAR  imaging;  they  can  be  expanded  in  Taylor 
series  as 

f  R(t)  =  R0  +  vt  +  \v't2  +  jv"f3  +  •  •  • 

l  6(0  =  Qf  +  ±fi'f2  +  \Wt3  +  •  •  • 

Taking  the  first  three  terms  of  the  Taylor  series,  we 
have 


A4 

f(x,o  =  XAexP 


k=[ 


(R0  +  (v  +  ykQ)t 


+  j(v'  +  xkQ2  +  yk£l')t2  +  ■■■)  . 

(3) 


The  simplified  phase  term  can  be  considered  as  a  sum 
of  polynomials  in  time  where  the  coefficients  of  each 
polynomial  are  determined  by  the  coordinates  of  the 
scatterers  and  the  target  motion  parameters  such  as 
velocity  and  angular  velocity.  The  constant  phase  term 
in  (3)  has  nothing  to  do  with  the  imaging  procedure 
and  thus  can  be  ignored.  While  the  first  order  term 
is  the  way  to  achieve  cross-range  resolution,  the 
second  and  higher  order  terms  are  the  source  of 
the  phase  error  which  cause  the  image  blurring. 

Hence  to  estimate  and  eliminate  these  higher  order 
phase  terms  is  the  main  task  of  the  following  motion 
compensation. 


III.  ADAPTIVE  JOINT  TIME-FREQUENCY  TECHNIQUE 

The  AJTF  technique  used  here  is  basically  a 
dechirping  procedure  just  like  other  conventional 
time-frequency  methods  for  detecting  a  chirp  signal 
in  the  noise  [15].  Nevertheless,  for  a  special  kind 
of  signal  such  as  radar  Doppler  frequency,  AJTF 
technique  proves  to  be  a  more  effective  and  efficient 
approach  to  track  the  Doppler  frequency  shifts  which 
behave  like  a  bunch  of  linear  or  higher  order  chirps  in 
the  time-frequency  plane.  For  simplicity,  we  illustrate 
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the  procedure  of  performing  AJTF  search  only  for 
linear  chirps.  As  shown  in  Fig.  1,  if  only  quadratic 
phase  error  terms  exist  in  (3),  the  radar  Doppler 
frequency  shifts  via  dwell  time  can  be  described  by 
those  linear  chirps  with  different  displacements  and 
slopes.  Each  chirp  represents  an  individual  point 
scattered  The  goal  of  AJTF  is  to  map  out  several 
strong  ones  and  extract  their  displacement  and  slope 
parameters.  For  this  goal,  a  set  of  basis  functions  are 
constructed  as. 


h  At)  =  exp[-j2n(f()t  +  i/,/2  +  •••)]■ 


(4) 


They  are  basically  a  collection  of  unit  chirps  with 
all  possible  displacement  and  slope  values,  (the 
dotted  line  in  Fig.  1).  Hence  projecting  the  radar 
Doppler  signal  to  these  unit  chirps  and  using  a  search 
procedure  to  maximize  the  projection  value,  the  basis 
function  which  looks  most  “alike”  to  the  strongest 
chirp  signal  will  be  picked  out.  The  projection 
procedure  is  searching  the  parameters  (/0,/j,...)  which 
satisfy 


\BP\2  =  max 
yoJi  — 


j  f(x,t)h‘p(t)dt\ 


(5) 


where  Bp  is  the  complex  coefficient  of  the  strongest 
chirp.  Therefore  we  can  extract  this  strongest 
chirp  signal  out  from  f(x,t)  and  running  the 
above  procedure  again  to  search  for  the  second 
strongest  chirp  in  the  signal.  For  the  specific  motion 
compensation  application  in  this  work,  we  only 
need  search  the  first  strongest  ccatterer  in  one  range 
cell,  then  look  for  other  scatterers  in  other  range 
cells  if  more  than  one  scatterer  needs  to  be  tracked. 
The  search  of  the  first-order  coefficient  can  be 
accomplished  by  using  fast  Fourier  transform  (FFT) 
algorithm,  then  only  a  one-dimensional  parameter 
search  is  required  to  find  /p  the  quadratic  phase  error. 
The  algorithm  can  also  be  extended  for  the  third  or 
higher  order  phase  errors  except  more  computation 
time  is  needed.  The  resolution  between  chirps  is 
expected  to  be  the  Fourier  resolution  of  the  full 
azimuth  aperture  resolution.  It  doesn’t  require  that  the 
tracked  point  scatterer  is  dominant  or  well  isolated  in 
the  range  cell. 
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Fig.  1.  Search  procedure  of  AJTF  processing  for  radar  signal 
with  quadratic  phase  error. 


motion  error  and  the  range  dependent  terms  represent 
wide-angle  imaging  error.  Thus  it  is  possible  to 
determine  all  the  motion  parameters  using  the 
information  provided  by  three  reference  scatterers. 
The  basic  idea  of  the  multiple  PPP  algorithm  is, 
suppose  that  there  are  three  prominent  point  scatterers 
whose  phase  components  are  known,  the  translational 
motion  error  can  be  removed  by  unwrapping  the 
higher  order  phase  component  of  the  first  prominent 
point.  Then  a  nonlinear  relationship  between  the 
rotation  angle  and  dwell  time  can  be  obtained  using 
the  phase  information  of  the  second  prominent 
points.  From  this  relationship,  the  uniformly  azimuth 
distributed  data  can  be  interpolated  from  the  original 
collected  data  in  nonuniform  azimuth  mesh  and 
thus  the  rotational  motion  error  can  be  removed.  At 
last,  the  rotation  rate  can  be  estimated  by  measuring 
the  phase  of  the  third  prominent  points,  so  that  the 
polar  reformatting  can  be  applied  to  eliminate  the 
wide-angle  imaging  error.  The  limitation  of  PPP  is 
assuming  the  availability  of  target  prominent  points 
so  that  they  can  be  easily  separated  out  from  the 
other  scatterers  using  a  windowing  technique.  The 
selection  and  tracking  of  those  points  is  often  done 
manually.  Here,  instead  of  an  interactive  procedure, 
AJTF  technique  is  used  as  a  powerful,  automatic  tool 
to  select  and  track  the  reference  points.  As  to  detailed 
procedure  to  conduct  PPP,  one  can  refer  to  [13,  14]. 


IV.  MODEL  FOR  MOTION  ERROR  ELIMINATION 

Once  the  phase  information  of  individual  scatterers 
is  measured.  The  next  task  is  to  eliminate  those 
quadratic  or  higher  order  phase  terms  for  all  the 
scatterers.  Here  we  implement  the  model  of  multiple 
PPP  algorithm  to  eliminate  the  phase  error  step  by 
step. 

Considering  the  right-hand  side  of  the  formula  (3), 
the  coefficients  of  quadratic  phase  terms  consist  of 
three  parts.  The  range  and  cross-range  independent 
terms  represent  translational  motion  error,  while  the 
cross-range  dependent  terms  represent  rotational 


V.  SOME  EXAMPLES 

We  use  the  simulated  step  frequency  data  from  a 
B727  airplane  as  an  example  to  describe  the  whole 
motion  compensation  procedure.  The  center  frequency 
of  radar  is  9  GHz  and  the  bandwidth  is  150  MHz. 

The  total  number  of  pulses  is  256.  After  compressing 
the  image  in  both  range  and  cross-range  direction, 
a  total  of  64  range  cells  and  256  cross-range  cells 
image  is  produced.  Fig.  2(a)  is  the  ISAR  image  after 
range  alignment  and  cross-range  centroid  focusing. 
Since  the  target  contains  a  significant  amount  of 
nonuniform  rotational  motion,  only  a  part  of  the 
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Fig.  2.  (a)  Simulated  B727  ISAR  image  after  range  alignment  and  coarse  motion  compensation,  (b)  STFT  spectrogram  of  signal  in 

particular  range  cell. 


airplane  is  well  focused  while  the  part  away  from  the 
focus  are  seriously  blurred.  To  describe  the  underlying 
mechanism  of  autofocusing  more  explicitly,  we  plot 
the  variation  of  Doppler  frequency  versus  dwell 
time  in  a  particular  range  cell  in  the  time-frequency 
plane  in  Fig.  2(b).  The  short  time  Fourier  transform 
(STFT)  is  used  to  generate  the  time-frequency 
spectrogram.  We  can  almost  see  a  bunch  of  lines 
in  the  spectrogram,  where  each  line  represents  the 
time-varying  phase  characteristics  of  a  scattering 
center  in  the  range  cell.  For  example,  the  displacement 
of  the  line  represents  the  coefficient  of  the  linear 
phase  terms  and  the  slope  of  the  line  represents  the 


coefficient  of  the  quadratic  phase  term.  Unfortunately, 
it  is  not  possible  to  extract  the  desired  parameters 
from  such  a  fuzzy  image  due  to  the  low  resolution 
of  the  STFT.  Instead,  the  AJTF  technique  is 
employed  because  of  its  full  azimuth  aperture 
resolution. 

First  we  search  the  first  reference  point  to 
be  focused  in  this  particular  range  cell.  AJTF 
automatically  pick  out  the  strongest  scatterer,  which 
is  represented  by  the  strongest  line  in  the  spectrogram. 
Multiplying  a  phase  correction  vector  to  the  original 
signal,  the  image  obtained  at  this  stage  is  shown  in 
Fig.  3(a).  As  seen  in  the  time-frequency  display  of 
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Fig.  3.  (a)  ISAR  image  obtained  by  focusing  on  first  reference  point,  (b)  STFT  spectrogram  of  signal  in  particular  range  cell. 


Fig.  3(b),  the  strongest  line  has  been  straightened 
and  shifted  to  the  center  of  the  cross-range  axis. 

This  step  is  only  for  translational  motion  removal 
and  thus  can  be  applied  directly  to  the  range  aligned 
data  in  place  of  conventional  translational  motion 
compensation  algorithm.  For  this  example,  it  is  used 
as  a  finer  tuning  to  the  previous  focusing  algorithm 
and  to  prepare  a  fixed  origin  for  the  next  steps.  As 
we  can  see  in  Fig.  3(a),  the  image  around  the  origin 
is  well  focused  while  the  other  part  is  still  defocused 
due  to  the  rotational  error.  Thus  a  second  reference 
point  must  be  selected  to  determine  the  rotation  rate 
nonuniformity  of  the  target.  The  second  reference 
point  is  selected  by  searching  the  strongest  scatterer 


in  another  range  cell  using  AJTF.  After  interpolating 
the  nonuniform  azimuth  mesh  to  uniform  azimuth 
mesh,  the  image  becomes  well  focused  everywhere 
on  the  target  as  in  Fig.  4(a).  As  we  can  see  in  the 
time-frequency  display  of  Fig.  4(b),  all  the  lines  are 
straightened,  implying  that  all  the  quadratic  phase 
terms  are  eliminated  and  all  the  points  are  well 
focused.  In  Fig.4  (c),  we  show  the  image  from  the 
simulated  data  without  any  added  motion  error  as 
a  reference  of  comparison.  For  this  example,  since 
two  point  focusing  has  made  the  image  well  focused, 
which  means  the  wide-angle  imaging  error  is  not 
serious,  there  is  no  necessity  to  do  the  third  point 
focusing. 
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Fig.  4.  (a)  Final  motion  compensated  B727  ISAR  image  by  focusing  on  two  reference  points,  (b)  STFT  spectrogram  of  signal  in 

particular  range  cell,  (c)  Simulated  B727  ISAR  image  without  any  motion  added  in. 


Another  example  is  a  simulated  MIG25  airplane. 
The  radar  center  frequency  is  9  GHz  and  bandwidth 
is  500  MHz.  A  total  of  64  range  cells  and  256 
cross-range  cells  are  used  for  the  imaging.  The 
coarsely  motion  compensated  image  is  in  Fig.  5(a). 
After  one  point  focusing  algorithm  like  above,  the 
focus  is  shifted  to  about  the  center  of  the  target,  but 
it  is  still  smeared  away  from  the  center,  as  Fig.  5(b). 
Consequently,  the  second  point  focusing  algorithm 
is  applied  but  doesn’t  improve  at  all,  which  means 
the  target  itself  doesn’t  contain  much  rotational  error. 
The  smearing  is  simply  because  the  large  rotation 
angle  makes  the  polar  reformatting  necessary.  Thus, 
a  third  reference  point  is  tracked  in  the  same  way  and 
its  phase  measured.  The  quadratic  phase  component 
provides  a  reasonable  estimate  of  the  angular 
rotation  rate.  For  this  target,  the  total  rotation  angle 


is  estimated  to  be  about  20  deg.  Once  the  angular 
spacing  is  found,  the  original  data  is  polar  reformatted 
to  eliminate  the  range  walk  and  cross-range  smearing 
caused  by  wide-angle  imaging.  The  final  image 
is  obtained  in  Fig.  5(c).  All  the  scatterers  are  well 
focused  in  their  individual  range  and  cross-range 
cells. 

Furthermore,  AJTF  technique  has  also  been 
utilized  for  motion  compensation  of  the  measurement 
data  from  a  few  real  airplanes.  It  proves  to  be  rather 
stable  and  insensitive  to  the  background  clutter. 


VI.  CONCLUSIONS 

An  AJTF  technique  based  on  time-frequency 
dechirping  concept  has  been  proposed  to  automatically 
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Fig.  5.  (a)  Simulated  MIG25  ISAR  image  before  motion  compensation,  (b)  ISAR  image  obtained  by  focusing  on  first  reference  point, 

(c)  Final  motion  compensated  ISAR  image  by  focusing  on  three  reference  points. 


perform  Doppler  tracking  to  single  or  multiple 
scatterers  in  ISAR  imaging  for  the  moving  target. 
After  applying  AJTF,  several  reference  point 
scatterers  are  automatically  selected  and  their 
phase  components  measured.  The  PPP  model 
thus  is  applied  for  the  compensation  of  both 
translational  and  rotational  motion  error,  as  well 
as  the  azimuth  spacing  estimation.  The  algorithm 
has  been  applied  to  two  simulated  airplanes  to 
illustrate  the  feasibility  and  the  detailed  procedure 
of  motion  compensation  and  automatically  polar 
reformatting.  The  high  resolution  and  stability, 
and  the  capability  to  track  individual  scatterers 
prove  to  be  the  main  advantages  of  AJTF  Doppler 
tracking. 
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A  Fast  Algorithm  for  Simulating  Doppler 
Spectra  of  Targets  with  Rotating  Parts  Using 
the  Shooting  and  Bouncing  Ray  Technique 

Raj  an  Bhalla  and  Hao  Ling 


Abstract — We  present  a  fast  algorithm  for  the  shooting  and  bouncing 
ray  technique  to  simulate  the  Doppler  spectra  of  targets  with  rotating 
parts.  The  fast  algorithm  gives  an  order  of  magnitude  speedup  in 
computation  time  while  maintaining  good  agreement  with  the  brute- force 
results. 

Index  Terms — Doppler  effect,  geometrical  theory  of  diffraction,  rotating 
bodies. 


I.  Introduction 

It  is  well  known  that  when  rotating  parts  exist  on  a  target  such 
as  propellers  or  rotors,  additional  features  are  introduced  in  the 
Doppler  frequency  spectra.  These  features  give  valuable  information 
about  the  structure  and  motion  of  the  rotating  part  and  can  be 
utilized  for  target  identification  applications  [l]-[3].  Therefor  d  is 
important  to  develop  the  capability  to  simulate  this  effect  from  first- 
principle  electromagnetics  and  to  relate  the  observed  effects  to  the 
target  structure  [4],  [5].  In  this  paper,  we  present  a  fast  method  to 
compute  the  Doppler  spectra  of  complex  targets  using  the  shooting 
and  bouncing  ray  (SBR)  technique  [6]-[8].  The  standard  methodology 
to  generate  the  Doppler  spectra  requires  the  computation  of  the 
scattered  field  at  different  time  snapshots  of  the  target  as  the  moving 
part  undergoes  motion.  The  resulting  time-varying  scattered  field  is 
then  Fourier  transformed  to  arrive  at  the  Doppler  spectrum.  From  the 
SBR  simulation  perspective,  each  time  snapshot  on  the  target  with 
a  different  moving  part  position  is  a  new  electromagnetic  problem 
and  requires  that  the  whole  computation  be  carried  out  from  scratch. 
This  can  be  very  computationally  intensive  for  a  complex  target.  For 
example,  a  typical  helicopter  rotor  having  a  blade  length  of  20  ft  and 
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otating  at  350  rpm  must  be  sampled  in  time  at  an  interval  of  0.34 
ns  or  less  to  avoid  aliasing  in  the  Doppler  spectrum.  In  order  to  map 
)ut  one  complete  rotor  revolution  at  10  GHz.  approximately  5110 
:ime  snapshots  of  the  target  are  needed.  If  we  use  SBR  to  carry  out 
he  scattering  calculation  at  each  rotor  position,  it  would  take  roughly 
100  days  (5110  x  30  min/run  on  a  typical  workstation)  to  compute 
he  Doppler  spectra  at  just  one  look  angle.  In  this  work,  we  present 
a  fast  Doppler  simulation  algorithm  for  the  SBR  method.  The  basic 
idea  is  to  take  advantage  of  the  instantaneous  Doppler  information 
available  for  each  ray  being  traced  through  the  target.  By  adjusting 
for  the  differential  phase  difference  due  to  target  motion  on  a  ray-by¬ 
ray  basis,  it  is  possible  to  generate  the  time  snapshots  of  the  target  at 
M  nearby  time  instances  from  ray  tracing  at  only  one  time  snapshot. 
Thus,  we  achieve  a  saving  of  M  ray-tracing  operations.  A  similar 
idea  has  been  devised  earlier  by  us  to  the  angular  extrapolation  of 
range  profiles  [8]. 


Fig.  1.  Each  time  a  ray  hits  a  moving  part,  it  undergoes  a  Doppler  frequency 
shift  depending  on  the  velocity  of  the  part,  vn  and  the  directions  of  the 
incident  and  reflected  rays  ktnt  and  kref. 


II.  Algorithm 

Based  on  SBR,  the  time-varying  scattered  field  of  a  target  contain¬ 
ing  moving  pans  takes  on  the  form 

E‘(t)  =  ]T  7 i(t)e~JViW  (1) 

t  r&yi 

where  7t  is  magnitude  contribution  of  each  ray  to  the  total  scattered 
field.  px  =  (27r/o/c)dt  is  the  phase  of  each  ray,  where  /o  is  the 
radar  frequency,  c  is  the  speed  of  light  and  dt  is  the  distance  traveled 
by  the  ray.  We  will  assume  that  the  target  is  stationary  except  for  the 
moving  parts.  As  time  evolves  the  moving  parts  undergo  motion  and 
those  rays  which  interact  with  the  moving  parts  will  have  a  different 
and  <pt  at  each  time  instance. 

We  shall  now  attempt  to  extrapolate  the  scattered  field  to  a  new 
time  t  from  information  at  to  by  making  several  approximations.  We 
first  assume  that  the  ray  hit  points  on  the  target  remain  the  same  at 
the  new  time  instance.  In  addition,  we  assume 


7»(f)  —  7«(fo) 


^(t)~cpl(f0)  + 


dip  t 
dt 


(t-to). 

*0 


(2a) 

(2b) 


Clearly,  these  approximations  hold  only  if  the  moving  parts  have  not 
undergone  large  displacements  from  to  to  t.  When  these  approxima¬ 
tions  are  valid,  however,  we  see  that  we  can  easily  compute  E*(t) 
from  the  available  ray  information  {7*  ( ),  p » ( to ) } ,  plus  the  quantity 
{dipt/dt^tQ  for  each  ray. 

To  compute  (d<pi/dt)  |to  for  each  ray,  it  is  more  convenient  to 
define 


foi  = 


1  dip  t 

2tt  dt 


*0 


(3) 


which  can  be  interpreted  as  the  “instantaneous”  Doppler  frequency 
of  ray  i.  As  shown  in  Fig.  1,  at  each  bounce  on  the  target,  a  ray 
undergoes  a  positive  or  negative  Doppler  shift  depending  on  the 
velocity  of  the  target  at  the  hit  point  v„  and  the  directions  of  the 
incident  and  reflected  rays  k;nc  and  kref.  If  we  assume  the  incident 
ray  is  at  a  frequency  the  frequency  of  the  reflected  ray  after 

ihe  impact  is  given  by 


/{")  _  1  _  5*1 


kjnc  *  Vri 


n  =  1.  2.  -V. 


jkref  •  Vn  ^ 


(4a) 


\fter  .V  bounces,  the  total  Doppler  frequency  shift  experienced  by 
he  ray  is  then 

'n,  =  f'y)  -  fiQ)  r4M 


where  =  fo  is  the  original  radar  frequency.  Therefore,  once  we 
have  computed  the  instantaneous  Doppler  experienced  by  each  ray 
via  (4),  we  have  the  necessary  information  to  carry  out  the  scattered 
field  extrapolation  to  other  time  instances. 

We  should  point  out  that  the  ray  summation  process  of  (1)  can  be 
accelerated  via  a  fast  Fourier  transform  (FFT)-based  fast  algorithm. 
The  details  of  the  algorithm  has  already  been  reported  in  [7]  and 
will  not  be  repeated  here.  The  result  of  applying  the  fast  algorithm 
is  that  the  ray  summation  takes  essentially  negligible  time  relative  to 
the  ray-tracing  time.  Therefore,  the  total  computational  time  speedup 
of  using  a  single  ray  trace  and  extrapolating  the  data  to  M  nearby 
time  instances  versus  generating  the  time  snapshots  via  brute  force 
is  the  factor  M. 


in.  Numerical  Results 

The  Doppler  spectrum  of  a  model  helicopter  with  a  rotor  having 
four  blades  is  generated  using  the  extrapolation  scheme  and  compared 
against  the  brute- force  approach.  The  target  is  a  model  helicopter.  The 
blade  length  of  the  rotor  is  10  ft.  One  commonly  used  display  of  the 
time-varying  Doppler  information  is  the  joint  time-frequency  spec¬ 
trogram.  It  is  generated  by  applying  a  short-time  Fourier  transform  to 
the  time-varying  scattered  field  data.  Fig.  2(a)  shows  the  joint  time- 
frequency  plot  of  the  helicopter  Doppler  spectrum  using  the  standard 
brute-force  approach.  For  this  example,  it  is  assumed  that  the  fuselage 
is  stationary  and  only  the  blades  are  rotating  at  350  rpm.  The  incident 
angle  on  the  target  is  10°  from  nose-on  and  the  radar  frequency  is 
2  GHz.  The  1400  time  samples  of  the  scattered  field  are  computed 
over  one  complete  revolution  of  the  rotor  blade  (360°  or  0.171  s). 
A  35-point  running  window  in  time  is  then  used  to  generate  the 
joint  time-frequency  plot.  The  total  computation  time  to  generate  the 
data  is  1400  x  2  min/run  =  2800  min  on  an  SGI  Indigo2  workstation 
(R4400/200MHz  CPU).  Fig.  2(b)  shows  the  joint  time-frequency  plot 
generated  using  the  extrapolation  scheme.  Ray  tracing  is  carried  out 
at  every  9°  of  rotor  revolution  and  the  rest  of  the  data  are  extrapolated 
from  the  ray  information  at  these  positions.  Therefore,  to  generate  the 
complete  joint  time-frequency  plot,  ray  tracing  is  carried  out  40  times. 
The  total  computation  time  for  the  fast  scheme  is  40  x  2  min/run  = 
80  min.  This  speedup  is  essentially  equal  to  the  theoretical  estimate 
of  35(=  1400/40)  since  the  extrapolation  overhead  is  negligible.  The 
joint  time-frequency  plot  generated  using  the  extrapolation  algorithm 
compares  well  with  that  from  the  brute-force  approach.  The  joint 
time-frequency  plots  show  four  fairly  distinct  patterns  corresponding 
to  the  Doppler  features  from  the  four  blades.  The  blade  flashes  at 
the  cardinal  angles  can  be  seen  for  each  blade.  As  the  blade  rotates, 
the  tin  mans  out  a  sinusoidal  nnttem  in  time  (or  aneleV  The  bodv 
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Fig.  2.  Comparison  of  the  joint  time-frequency  plots  generated  using 
the  standard  brute-force  approach  and  the  fast  extrapolation  algorithm,  (a) 
Brute-force  approach,  (b)  Fast  extrapolation  scheme. 

\ 

of  the  helicopter  corresponds  to  the  zero  Doppler  bin  since  it  is 
assumed  to  be  stationary.  There  is  some  spreading  of  scattered  energy 
about  the  zero  Doppler  bin  due  to  body-rotor  interactions.  Finally,  we 
should  point  out  that  the  9°  extrapolation  window  used  in  the  present 
example  is  chosen  by  experience  and  will  in  general  be  governed  by 
target  complexity. 
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ceivcd  signal.  For  a  moving  target,  the  mea¬ 
surement  of' the  target's  velocity  is  based  on  the 
well-known  Doppler  effect.  If  the  ra¬ 
dar-transmitted  signal  is  at  frequency/;  the  re¬ 
flected  signal  from  the  moving  target  is 
subjected  to  a  Doppler  frequency  shift  from  its 
transmitted  frequency,  /*+  //  This  frequency 
shift  is  induced  bv  relative  motion  between  the 
radar  and  the  target  [  2  ].  In  the  case  where  a  tar¬ 
get  has  a  radial  velocity  i\  the  Doppler  fre¬ 
quency  shift  is  determined  bv  the  radial 
velocity  of  the  target  and  the  wavelength  of  the 
radar  transmitted  waveform 


A  /.  Example  of  radar  operational  scenario:  transmitter ,  receiver,  target,  noise, 
and  clutter. 


fd  -~2f- 


(3) 


▲  2.  Examples  of  manmade  targets:  surface,  edge,  corner, 
trihedral,  and  cavity. 


where  R  is  the  range  from  the  radar  to  the  point  target 
and  c  is  the  wave  propagation  velocity. 

When  the  target  moves  with  a  velocity  v  relative  to  the 
radar,  called  the  radial  velocity,  the  radar  signal  must 
travel  an  extra  longer  or  shorter  distance  to  reach  the  tar¬ 
get.  Therefore,  the  signal  received  at  time  t  is  reflected 
from  the  target  at  time  (f-T  (t)/2),  and  the  received  signal 
becomes  sR(t)(jpsr[t-x(t)U  where  the  round-trip  travel 
time  is  a  time- varying  delay  [2], [7]. 

In  addition  to  the  signal  reflected  bv  the  target,  there  is 
also  additive  noise.  The  signal-to-noise  ratio  (SNR)  at  the 
radar  receiver  is  determined  by  the  received  power  re¬ 
flected  from  the  target,  the  noise  figure,  and  the  band¬ 
width  of  the  receiver.  Improvement  in  the  SNR  will 
increase  the  probability  of  target  detection  and  the  accu- 
raev  of  parameter  estimation.  A  radar  usually  transmits  a 
sequence  of  pulses  or  other  signal  waveforms  with  a  pulse 
repetition  frequency  (PRF)  usually  required  by  the  maxi¬ 
mum  range  of  detection.  Target  information  may  be  ex¬ 
amined  directly  from  the  radar-range  profile,  i.e.,  the 
distribution  of  target  reflectivity  along  the  radar  line  of 
sight  to  the  target,  or  from  its  frequency  spectrum  by  ap¬ 
plying  the  Fourier  transform  [8]-[ll]. 

In  general,  the  radar  processes  the  received  signal  and 
extracts  information  about  the  target.  The  range  to  the 
target,  i.e..  the  distance  from  the  radar  to  the  target  mea¬ 
sured  along  the  radar  line  of  sight,  is  estimated  by  the 
time-delav  between  the  transmitted  signal  and  the  re- 


Thus,  if  a  target  is  closing  on  the  radar  at  a  velocity  v  = 
-1000  (ft/s)  =  -304.8  (m/s),  the  Doppler  frequency  shift 
for  X-band  radar  at  9842  MHz  is  20  kHz. 

Radar  targets,  especially  man-made  targets,  usually 
can  be  considered  as  a  collection  of  point  scatterers.  These 
scatterers  may  have  a  wide  variety  of  reflecting  or  back- 
scattering  behaviors  [12],  [13].  They  can  include  sur¬ 
faces,  edges,  corners,  dihedrals,  trihedrals,  and  cavities 
(Fig.  2).  Each  type  of scatterer  has  a  different  backscatter- 
ing  behavior,  which  can  provide  a  wav  to  recognize  the 
target.  For  example,  when  cavities  or  duct-tvpe  structures 
are  present  in  targets,  multiple  bounces  appear  in  the  ra¬ 
dar  image  as  “clouds'1  extended  in  the  range  domain  [43] . 
However,  the  same  mechanism  mav  provide  important 
features  that,  if  properly  interpreted,  can  be  important 
factors  in  the  target  recognition  process.  In  the  next  sec¬ 
tion,  we  will  further  describe  radar  backscattering  and  in¬ 
troduce  time-frequenev  analysis  applications  to  radar 
backscattering  and  feature  extraction. 

Range  profiles  have  been  examined  for  target  recogni¬ 
tion,  especially  in  cases  when  a  radar  image  cannot  be 
generated  due  to  various  reasons,  such  as  short  radar 
dwell  time  on  the  target  [9].  With  high  range  resolution, 
the  range  profile  can  provide  target  information  about 
target  length  and  positions  of  stronger  scatterers,  such  as 
a  radar  dish,  engine  intakes,  and  other  strong  scattering 
centers  (Fig.  3).  The  dimension  transverse  to  the  radar 
line  of  sight  is  called  the  cross  range.  Because  the  Doppler 
shift  of  a  scatterer  in  a  target  is  proportional  to  the  cross 


A  3 .  Example  of  an  aircraft's  range  profile. 


range  of  the  scattercr  bv  a  scale  factor,  the  projection  of 
the  target  reflectivity  distribution  on  the  cross-range  di¬ 
mension  can  be  obtained  from  the  distribution  of  Dopp¬ 
ler  shifts,  called  the  Doppler  profile.  With 
high-resolution  Doppler  profiles,  the  locations  of  stron¬ 
ger  scatterers  and  the  target's  extents  in  the  Doppler  di¬ 
mension  can  be  obtained  as  illustrated  in  Fig.  4.  By 
combining  range  profiles  and  Doppler  profiles,  a 
2-D-radar  image  mav  be  generated  [ 1 3 ] - [  16].  The  radar 
image  is  a  spatial  distribution  of  the  target’s  reflectivity 
mapped  onto  a  range  and  Doppler  plane.  The 
range- Doppler  image  can  be  converted  to  this  2-D  image 
if  we  have  accurate  knowledge  of  the  scale  factor.  This  is 
determined  bv  the  rotation  rate,  the  distance  of  the  scat- 
terer  from  its  rotation  center,  and  the  wavelength  of  the 
transmitted  signal  [15],  [16]. 

An  important  factor  of  the  image  quality  is  its  resolu¬ 
tion.  Resolution  is  the  ability  to  separate  closely  spaced 
scatterers  in  range  and  in  cross  range.  The  minimum  dis¬ 
tance  in  the  range  Ar  combined  with  in  the  cross  range  or 
azimuth  A.;.  bv  which  two  point-scatterers  may  be  sepa¬ 
rated  is  the  resolution  of  the  image.  A  rectangle  with  sides 
A.  and  A is  called  a  resolution  cell.  Range  resolution  is  de¬ 
termined  bv  the  frequency  bandwidth  passed  by  the 
transmitter  and  receiver.  For  example,  a  500-MFIz  band¬ 
width  can  vield  a  range  resolution  of  1  ft.  The  width  of  the 
bandwidth  obtained  depends  on  the  radar  operating  fre¬ 
quency.  For  an  X-band  radar  operated  at  a  10,000-MHz 
frequency,  a  bandwidth  of  5%  of  the  radar  operation  fre¬ 
quency,  i.e.,  500  MHz,  is  possible.  Generally,  higher  res¬ 
olution  is  more  easily  achieved  in  the  range  domain  than 
in  the  cross-range  domain.  To  obtain  high  cross-range 
resolution  without  utilizing  an  unpractically  large  an¬ 
tenna  aperture,  synthetic  aperture  processing  is  widely 
used.  It  coherently  combines  signals  obtained  from  se¬ 
quences  of  small  apertures  at  different  aspect  angles  of  the 
target  to  emulate  a  result  which  could  be  obtained  from  a 
large  antenna  aperture  [14],  [15]. 

Coherent  processing  maintains  the  relative  phases  of 
successive  pulses.  Thus,  the  phase  from  pulse  to  pulse  is 
preserved  and  a  phase  correction  can  be  applied  to  the  re- 
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\  4.  An  illustration  of  the  Doppler  profile  of  an  aircraft . 


To  use  the  Fourier  transform 
properly,  it  is  assumed  that  the 
frequency  contents  of  the 
analyzed  signal  must  be 
time-invariant. 

turns  to  make  them  coherent  for  successive  interpulse  pe¬ 
riods.  If  radar  returns  are  processed  coherently,  the  pro¬ 
cessed  data  retains  both  the  amplitude  and  the  phase 
information  about  the  target.  The  amplitude  is  related  to 
the  radar  cross  section  of  the  target  { a  measure  of  the  abil¬ 
ity  to  reflect  electromagnetic  waves  >  and  the  phase  is  re¬ 
lated  to  the  radial  velocity  of  the  target.  The  time  history 
series  of  samples  taken  at  each  range  forms  the  basis  for 
radar  signal  and  image  processing. 

Imaging  of  moving  targets*  using  radar  has  been  a  ma¬ 
jor  challenge.  Techniques  for  high- resolution  radar  imag¬ 
ing  are  based  on  synthetic  aperture  formation.  When 
there  is  a  relative  motion  between  the  radar  and  the  tar¬ 
get,  a  synthetic  aperture  can  be  formed.  Traditionally, 
synthetic  aperture  radar  (SAR)  refers  to  the  situation 
w'here  the  radar  is  moving  [I4]-[16],  w'hile  inverse  syn¬ 
thetic  aperture  radar  (ISAR)  refers  to  the  situation  where 
the  target  is  moving  [13],  [15]. 

To  form  an  image  from  the  radar  data,  accurate  mo- 
Lion  compensation  is  needed.  For  ISAR,  if  the 
translational  motion  is  compensated  for  by  using  a  refer¬ 
ence  point  on  the  target  such  as  the  center  of  rotation, 
then  the  moving  target  becomes  a  rotating  target  around 
that  center,  and  forms  an  equivalently  larger  circular  aper¬ 
ture  focused  at  the  center.  The  target  can  be  considered  to 
be  a  set  of  individual  point  scatterers,  each  w  ith  a  radial 
velocity  or  Doppler  frequency  shift  to  the  radar.  Thus,  the 
distribution  of  the  radar  reflectivity  of  the  target  can  be 
measured  bv  the  Doppler  frequency  spectrum  at  each 
range  estimated  bv  taking  the  Fourier  transform  over  the 
observation  time  interval. 

To  use  the  Fourier  transform  properly,  it  is  assumed 
that  the  frequency  contents  of  the  analyzed  signal  must  be 
time-invariant.  With  this  assumption,  a  long  observation 
time  results  in  high  Doppler  resolution.  However,  when 
rhe  target  moves,  Doppler  frequency  shifts  are 
rime-varying,  and  the  assumption  of  time-invariant 
Doppler  frequency  shifts  is  no  longer  valid.  Thus,  the 
Doppler  spectrum  becomes  smeared,  degrading  the 
cross-range  resolution,  and  the  radar  image  becomes 
blurred.  There  are  main’  methods,  some  called 
autorocusing  and  others  called  motion-compensation, 
used  in  solving  the  Doppler  smearing  and  image  blurring 
problem  j  1  '”1-126  | .  Most  methods  are  Fourier-based  ap¬ 
proaches  that  attempt  to  flatten  Doppler  spectra  or  indi¬ 
vidual  scatterers  bv  using  sophisticated  preprocessing 
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A  5.  Electromagnetic  mechanisms  are  manifested  in  the  joint 
time-frequency  image  as  distinct  features,  (a)  Scattering  center 
(b)  resonance  (c)  material  dispersion,  and  (d)  structural  disper- 


approaches.  Others  use  modern  spectral  analysis  to 
achieve  sharper  images  with  shorter  data  samples 
[27]-[30].  Sets  of  radar  data  are  available  from  http://air- 
borne.nrLnavy.mil/~vchen/tftsa.html.  In  another  sec¬ 
tion,  we  introduce  an  image  formation  based  on 
time-frequency  analysis,  which  can  resolve  the  im¬ 
age-blurring  problem  without  resorting  to  sophisticated 
preprocessing  algorithms. 

Applications  of  JTF  Analysis  to 
Radar  Backscattering 

Time-Frequency  Analysis  of  Radar  Range  Profiles 

With  the  radar  operating  at  a  sufficiendv  high  bandwidth 
(i.e.,  equivalent  to  a  short  pulse),  the  resulting 
high-resoludon  range  profile  can  be  interpreted  as  a  map¬ 
ping  of  the  reflectivity  of  the  target  along  the  radar  line  of 
sight.  With  simple  targets,  a  range  profile  tvpicallv  con¬ 
sists  of  a  series  of  distinct  peaks  that  can  be  related  spa¬ 
tially  to  the  scattering  centers  on  the  target.  These 
features  are  often  examined  for  signature  diagnostic  and 
target  recognition  applications.  For  real  targets,  however, 
the  scattering  characteristics  are  usuallv  more  complex. 
For  example,  the  scattering  from  some  components  on  a 
target  can  be  strongly  dispersive  as  a  funcdon  of  fre¬ 
quency,  and  may  give  rise  to  extended  returns  in  the  range 
domain.  These  dispersive  scattering  phenomena  can  be 
difficult  to  interpret  from  the  time-domain  range  profile. 

JTF  methods  described  in  [1]  have  been  used  to  ana¬ 
lyze  electromagnetic  backscattered  data  with  good  suc¬ 
cess  [31] -[42].  The  JTF  representation  of  a  signal  is  a 
2-D-feature  space  that  facilitates  the  visualization  and  in¬ 
terpretation  of  complex  electromagnetic  wave  phenom¬ 
enology.  In  this  feature  space,  discrete  time  events  such  as 
scattering  centers,  discrete  frequency  events  such  as  reso¬ 
nances,  and  dispersive  mechanisms  due  to  surface  waves 
and  guided  modes  can  all  be  displaved  simultaneousiv. 
This  can  lead  to  more  insight  into  the  complex  electro- 


A  6.  Joint  time-frequency  (JTF)  image  of  the  backscattering  data 
from  a  coated  plate  with  a  gap  in  the  coating.  The  JTF  image  is 
generated  by  the  short-time  Fourier  transform.  Those  features 
which  show  slanting  in  the  JTF  plane  are  associated  with  the 
dispersive  surface  wave  mechanisms  in  the  coating. 

magnetic  wave  propagation  and  scattering  mechanisms 
than  the  information  available  from  either  time  or  fre¬ 
quency  domain  representations  alone. 

Shown  in  Fig.  5  are  the  time-frequency  features  of 
some  commonly  encountered  scattering  mechanisms.  A 
discrete  event  in  time,  for  example,  could  be  due  to  wave 
scattering  from  a  spatially  localized  scattering  center  on  a 
structure.  It  shows  up  as  a  vertical  line  [Fig.  5(a)]  in  the 
image  because  it  occurs  at  a  particular  instance,  but  over 
all  frequencies.  A  target  resonance,  e.g.,  the  return  from  a 
partially  open  cavity,  is  a  scattering  event  that  becomes 
prominent  at  a  particular  frequency.  It  shows  up  as  a  hori¬ 
zontal  line  in  the  JTF  plane  [Fig.  5(b)].  Dispersive  phe¬ 
nomena,  on  the  other  hand,  are  characterized  by  slanted 
curves  in  the  time-frequency  image.  For  instance,  surface 
wave  mechanisms  due  to  material  coatings  are  character¬ 
ized  by  curves  with  a  positive  slope  [Fig.  5(c)].  Another 
type  of  dispersion  arises  from  waveguide  structures. 
These  “structural  dispersion”  mechanisms  are  character¬ 
ized  by  curves  with  a  negative  slope  in  the  time-frequency 
image  [Fig.  5  (d)  ].  All  of  the  above  phenomena  have  been 
observed  in  a  wide  variety  of  structures,  from  simulation 
data  on  canonical  structures  to  measurement  data  on 
complex  platforms  [31]-[42].  Below,  two  examples  are 
presented  to  demonstrate  the  unique  features  of  electro¬ 
magnetic  scattering  mechanisms  in  the  JTF  plane. 

In  the  first  example,  the  backscattered  data  from  a  di¬ 
electric-coated  plate  with  a  gap  in  the  coating  is  consid¬ 
ered  [35]  (Fig.  6).  The  radar  signal  is  incident  edge-on 
from  the  left,  with  the  incident  electric  field  polarized  in 
the  vertical  direction.  The  radar  frequency  response  has 
been  generated  by  computer  simulation  and  verified  by 
laboratory  measurement.  The  simulation  result  is  shown 
along  the  vertical  frequency  axis.  The  time-domain  re¬ 
sponse,  or  equivalently,  the  range  profile  (where  range  = 
ct/2),  is  obtained  by  Fourier  transforming  the 
1 .7-to-18-GHz  data.  The  resulting  range  profile  is  shown 
along  the  horizontal  time  axis.  It  appears  that  three  dis¬ 
tinct  pulses  are  present.  However,  the  second  and  third 
pulses  are  spread  out  in  range. 
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In  order  to  resolve  in  finer  detail  the  dispersive  scatter¬ 
ing  mechanisms  in  this  coated  plate,  the  spectrogram  of 
the  backscattered  signal  is  generated  using  the  short-time 
Fourier  transform.  As  can  be  seen,  the  scattering  mecha¬ 
nisms  are  much  more  apparent  in  the  2-D  JTF  plane  than 
in  either  the  time  or  the  frequency  domain.  In  particular, 
it  is  observed  that  the  third  broad  pulse  in  the  time  do¬ 
main  actually  consists  of  three  separate  scattering  mecha¬ 
nisms  (labeled  as  3a,  3b,  and  3c).  As  the  frequency 
approaches  zero,  the  propagation  delays  of  mechanisms 
3a,  3b,  and  3c  approach  the  same  value.  As  the  frequency 
increases,  the  pulses  have  different  propagation  delays, 
and  become  clearly  separated. 

Slanted  curves  in  the  time-frequency  plane  (like  mech¬ 
anisms  2,  3b,  and  3c)  are  characteristic  of  dispersive  be¬ 
havior.  In  the  case  of  the  coated  plate,  surface  waves 
excited  in  the  coating  give  rise  to  the  dispersive  mecha¬ 
nisms.  At  frequencies  well  above  cutoff,  the  surface  wave 
is  tightlv  bound  to  the  dielectric  and  the  wave  velocity  ap¬ 
proaches  the  slow  dielectric  velocity.  Near  cutoff,  the  sur¬ 
face  wave  velocity  approaches  that  of  free  space  and 
exhibits  a  shorter  propagation  delay.  Therefore,  in  the 
time-frequency  plane,  the  surface  wave  phenomena  show 
up  as  slanted  curves  with  a  positive  slope. 

Based  on  the  propagation  delay  considerations  and  the 
above  observation,  it  is  possible  to  pinpoint  the  five  domi¬ 
nant  scattering  mechanisms.  They  are  shown  on  the  left  in 
Fig.  6,  which  clearly  indicates  that  mechanisms  2,  3b,  and 
3c  include  surface  wave  propagation.  For  the  polarization 
under  consideration  and  the  frequency  range  of  the  data, 
the  TMq  surface  wave  mode  that  has  zero  cutoff  is  the  only 


mode  that  can  propagate  in  the  dielectric  (the  TMj  mode 
has  a  cutoff  frequency  of  23.4  GHz) .  Finally,  higher-order 
scattering  mechanisms  can  be  observed  during  the 
late-time  portion  of  the  data,  but  they  are  very  weak. 


A  discrete  event  in  time,  for 
example,  could  be  due  to  wave 
scattering  from  a  spatially 
localized  scattering  center 
on  a  structure. 

In  the  second  example,  the  scattering  from  a  slotted 
waveguide  structure  is  considered  [39].  The  geometry  is 
shown  in  Fig.  8,  where  a  long  rectangular  waveguide  is 
flush  mounted  in  a  conducting  ground  plane.  Two  nar¬ 
row  slots  are  opened  on  each  end  of  the  ground  plane. 
The  structure  is  excited  by  a  horizontally  polarized  radar 
signal  occurring  at  an  angle  of  30°  with  respect  to  the  ver¬ 
tical.  The  backscattered  data  are  generated  by  a  computa¬ 
tional  electromagnetic  simulator  based  on  the  method  of 
moments.  The  data  are  generated  from  0.025  to  10  GHz 
in  25-MHz  increments. 

The  time-frequency  representation  of  the  data,  ob¬ 
tained  using  the  short-time  Fourier  transform,  is  shown 
on  the  right  in  Fig.  7.  The  two  early-time  vertical  lines 
correspond  with  the  exterior  scattering  centers  from  the 
slots.  The  other  curves  are  related  to  signals  coupled  into 


▲  7.  Joint  time-frequency  image  of  the  backscattering  data  from  a  slotted  waveguide  structure.  The  data  were  simulated  using  a 
method  of  moments  solver.  The  JTF  image  is  generated  by  the  short-time  Fourier  transform.  The  JTF  image  shows  both  the 
early-time  discrete-time  returns  from  the  slot  exterior  and  the  late-time  dispersive  mechanisms  due  to  modal  propagation 


inside  the  waveguide. 
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A  8.  Joint  time-frequency  processing  is  appiied  to  the  range  di¬ 
mension  of  the  conventional  range/cross-range  ISAR  image  to 
gain  an  additional  frequency  dimension.  By  examining  how 
the  resulting  images  vary  as  a  function  of  frequency,  the  fre¬ 
quency-independent  features  can  be  separated  from  the  fre¬ 
quency  dependent  ones  and  displayed  in  an  appropriate 
feature  space. 


the  waveguide  (Fig.  7,  again).  When  the  wave  reaches  the 
first  slot,  some  energy  is  coupled  into  the  waveguide, 
propagadng  to  the  other  end  as  a  sum  of  waveguide 
modes.  The  energy  carried  by  these  modes  begins  to 
reradiate  through  the  other  slot  after  a  time  given  by  L/c 
where  L  is  the  length  between  the  two  slots.  However, 
this  is  the  time  delay  only  for  frequencies  well  above  the 
modal  cutoff  for  which  the  modal  group  velocity  ap¬ 
proaches  c.  For  frequencies  approaching  the  cutoff  fre¬ 
quency  of  the  respective  mode,  the  group  velocity  tends 
toward  zero  and  the  time  delay  goes  to  infinity.  Conse¬ 
quently,  each  modal  dispersion  behavior  is  manifested  as  a 
time-frequency  trajectory  with  negative  slope  [as  illus¬ 
trated  earlier  in  Fig.  5(d)].  This  behavior  can  be  clearly 
identified  in  the  spectrogram,  where  the  presence  of  two 
modal  dispersion  curves  with  cutoffs  at  3  and  6  GHz  are 
observed.  They  correspond  to  the  TE10  and  TE20  modes  in 
the  waveguide.  Note  that  the  amplitude  variation  of  the 
signal  along  these  curves  is  governed  by  the  coupling 
mechanisms  through  the  slot  apertures  and  is  considerably 
more  complex. 

Because  multiple  reflections  occur,  we  also  see  other  dis¬ 
persion  curves  with  greater  delays  during  the  late-time  por¬ 
tions  of  the  return.  The  first  one  corresponds  with  the  energy 
that,  upon  reaching  the  other  end,  reflects  back  and  radiates 
through  the  slot  on  the  left.  The  next  is  the  three-bounce 
mechanism.  Note  that  energy  is  also  coupled  into  the  wave¬ 
guide  through  the  slot  on  the  right,  and  through  similar 
mechanisms,  generates  dispersion  curves  in  the 
time-ftequencv  plot  depending  on  the  number  of  bounces. 

From  the  above  two  examples,  it  can  be  seen  that  the 
JTF  representation  can  aid  in  the  interpretation  of  com¬ 


plex  electromagnetic  phenomena.  Furthermore,  the  JTF 
features  can  be  well  understood  in  terms  of  the  tar¬ 
get-scattering  phvsics.  For  radar  applications,  the 
time-frequency  representation  is  particularly  effective  in 
identifying  scattering  mechanisms  in  targets  containing 
“sub-skinline”  structures  such  as  inlet  ducts,  antenna  win¬ 
dows,  and  material  coatings. 


Physics-Based  Feature  Extraction  from 
Radar  imagery 

The  JTF  processing  of  the  one-dimensional  range  profile 
described  above  can  be  further  extended  to  deal  with  2-D 
radar  imaging.  ISAR  imaging,  as  discussed  earlier,  is  a  ro¬ 
bust  process  for  mapping  the  position  and  magnitude  of 
the  point  scatterers  on  a  target  from  multi- frequency  and 
multi-aspect  backscattered  data.  However,  for  complex 
targets  containing  other  scattering  phenomena  such  as 
resonances  and  dispersive  mechanisms,  image  artifacts 
are  often  encountered  in  the  resulting  ISAR  image.  One 
important  example  is  the  scattering  from  the  engine  in¬ 
let/exhaust  duct  on  aircraft.  It  is  a  dominant  contributor 
to  the  overall  scattering  from  the  target,  yet  its  wave¬ 
guide-like  structure  and  the  associated  fre¬ 
quency-dependent  scattering  mechanisms  make  it  a 
nonpoint-scattering  feature. 

When  processed  and  displayed  by  the  conventional 
ISAR  algorithm,  the  inlet  return  results  in  an  image  fea¬ 
ture  which  is  not  well-focused,  is  not  related  to  the  spatial 
location  of  the  scatterer,  and  can  often  obscure  other  im¬ 
portant  point  features  on  the  target.  Therefore,  it  would 
be  useful  to  automatically  remove  these  artifacts  from  the 
ISAR  image,  leading  to  a  clean  ISAR  image  containing 
only  physically  meaningful  point  scatterers.  Further¬ 
more,  the  inlet  features  extracted  can  be  better  displayed 
in  a  more  meaningful  feature  space  to  identify  target  reso¬ 
nances  and  cutoff  phenomena. 

JTF  processing  can  be  applied  to  ISAR  image  process¬ 
ing  to  accomplish  the  above  objective  [43].  The  concep¬ 
tual  idea  behind  the  JTF  ISAR  algorithm  is  to  apply  JTF 
processing  to  the  range  (or  time)  axis  of  the  conventional 
range/cross  range  ISAR  image  to  gain  an  additional  fre- 
quencv  dimension.  The  result  is  a  3D  range/ 
cross-range/frequency  matrix,  with  each 
range/cross-range  slice  of  this  matrix  representing  an 
ISAR  image  at  a  particular  frequency  (Fig.  8).  Conse¬ 
quently,  by  examining  how  the  ISAR  image  varies  with 
frequencv,  we  can  distinguish  the  frequency-independent 
scattering  mechanisms  from  the  frequency-dependent 
ones.  In  the  actual  implementation  of  the  JTF  ISAR,  the 
choice  of  the  JTF  processing  engine  is  critical  to  preserv¬ 
ing  range  resolution.  This  is  demonstrated  below  using 
the  adaptive  Gaussian  representation  proposed  by  Qian 
and  Chen  [44]. 

The  signal-adaptive  Gaussian  representation  is  a 
high- resolution  time-frequency  representation.  (A  simi¬ 
lar  algorithm  called  matching  pursuit  was  developed  in¬ 
dependently  at  around  the  same  time  by  Mallat  and 
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Zhang  [45]).  The  objective  of  this  method  is  to  expand  a 
signal  f(t)  in  terms  of  Gaussian-moduiated  exponential 
basis  functions  hp(t)  with  an  adjustable  standard  deviation 
ap  and  a  time-frequency  center  (tpfp): 


m^Bphp(t) 

p= i 


where 


be(t)=(na])^25cx p 


'  2°2>  J 


exp(j2nf  t). 


(4) 


(5) 


Note  that  the  basis  function  has  a  dual  form  in  its  Fou¬ 
rier  transform  representation: 
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Therefore,  these  basis  functions  have  a  time-frequency  ex¬ 
tent  given  byc;  and  (l/2ap),  respectively.  The  coefficients 
Bp  are  found  one  at  a  time  by  an  iterative  procedure.  One 
begins  at  the  stagey  =  1  and  chooses  the  parameters  ap  tp 
and  fp,  such  that  hp  (t)  is  most  “similar”  to  f(t),  that  is: 


NX™*  |f fH(t)b-p(t)dt\2 


(7) 

where f0(t)  =f(t).  For p  >  1  ,fp  (t)  is  the  remainder  after 
the  orthogonal  projection  of  fp_,  (, t)  onto  hp  (t)  has  been  re¬ 
moved  from  the  signal: 
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(c)  Enhanced  ISAR  Image 
After  JTF  Processing 

Enhanced  ISAR 
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(d)  Extracted  Inlet  Return  in 
the  Frequency  Aspect  Plane 
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9.  Adaptive  JTF  ISAR  processing  using  the  test  range  data  from  a  1:30  scaled  model  of  the  VFY-218  airplane,  (a)  Geometry  of  the 
VFY-218  with  two  deep  inlet  ducts,  (b)  Conventional  ISAR  image  for  frequency  from  8  to  16  GHz,  aspect  from  0°  to  40°.  The  dispersive 
cloud  over  the  right  wing  is  due  to  the  return  from  the  left  inlet  duct,  (c)  Enhanced  ISAR  image  obtained  by  removing  non-point  scat¬ 
tering  mechanisms  using  the  AJTFISAR  algorithm.  The  inlet  return  has  been  removed  from  the  image,  (d)  Frequency-aspect  display 
showing  the  extracted  features  from  the  engine  return.  Prominent  resonant  frequencies  are  observed  in  this  feature  space. 


MARCH  1999 


IEEE  SIGNAL  PROCESSING  MAGA7INF 


By  replacing  the  Fourier  trans¬ 
form  with  a  high-resolution 
time-frequency  transform,  the 
image  blurring  caused  by  the 
time-varying  Doppler  frequency 
shifts  can  be  mitigated  without 
applying  sophisticated  focusing 
algorithms. 

/,(*)=/,.,  (8) 

This  procedure  is  repeated  to  generate  as  many  coefficients 
as  needed  to  accurately  represent  the  original  signal. 

The  adaptive  Gaussian  representation  has  two  distinct 
advantages  over  conventional  time-frequency  techniques 
such  as  the  short-time  Fourier  transform.  First,  it  is  a 
parametric  procedure  that  results  in  very  high 
time-frequency  resolution.  More  importandv  for  our  ap¬ 
plication,  the  adaptive  spectrogram  allows  us  to  distin¬ 
guish  the  frequency-dependent  events  from  the 
frequency-independent  ones  automatically  through  the 
extent  of  the  basis  functions.  From  (5),  it  can  be  seen  that 
scattering  centers,  i.e.,  signals  witn  very  narrow  lengths  in 
time,  will  be  well  represented  by  basis  functions  with  very 
small  <5p.  Frequency  resonances,  on  the  other  hand,  will  be 
better  depicted  by  large  ap.  Therefore,  if  we  reconstruct 
the  ISAR  image  using  only  those  Gaussian  bases  with 
small  variances,  a  much  cleaner  image  can  be  obtained 
showing  only  the  scattering  centers.  The  remaining 
mechanisms,  i.e.,  those  related  to  the  large  variance 
Gaussians,  are  more  meaningful  when  viewed  in  a  dual 
frequencv-aspect  display,  where  resonances  and  other  fre¬ 
quency-dependent  mechanisms  can  be  better  identified. 

The  algorithm  is  demonstrated  by  using  the  chamber 
measurement  data  of  a  1:30  scale  model  Lockheed 
VFY-218  airplane  provided  by  the  Electromagnetic 
Code  Consortium  [46].  The  airplane  has  two  long  engine 
inlet  ducts,  which  are  rectangular  at  the  open  ends,  but 
merge  together  into  one  circular  section  before  reaching  a 
single-compressor  face  [Fig.  9(a)].  As  we  can  clearly  see 
in  the  conventional  ISAR  image  for  the  horizontal  polar¬ 
ization  at  20°  near  nose-on,  the  large  cloud  outside  of  the 
airframe  structure  is  the  inlet  return  [Fig.  9(b)].  Fig.  9(c) 
shows  the  enhanced  ISAR  image  of  Fig.  9(b),  obtained 
bv  applying  the  JTF  ISAR  algorithm  and  keeping  only 
the  small  variance  Gaussians.  We  see  that  only  the  scatter¬ 
ing  center  part  of  the  original  signal  remains  in  the  image 
as  expected.  Notice  that  the  strong  return  due  to  engine 
inlet  has  been  removed,  but  the  scattering  from  the  right 
wing  tip  remains.  Fig.  9(d)  shows  the  frequency-aspect 


display  of  the  high  variance  Gaussians.  A  number  of 
equispaced  vertical  lines  can  be  seen  between  10.5  and 
13.5  GHz.  Given  the  dimensions  of  the  rectangular  in¬ 
let  opening,  we  estimate  that  these  frequencies  corre¬ 
spond  approximately  to  the  second  cutoff  frequency  of 
the  waveguide-like  inlet.  This  information  may  be 
unique  to  the  particular  inlet  structure  under  consider¬ 
ation  and  mav  be  useful  as  an  additional  feature  vector 
in  target  classification  [47].  The  processed  images  for 
observation  angles  varying  from  0-360°  have  been 
made  into  an  MPEG  movie  that  can  be  downloaded 
from  http://lingO.ece.utexas.edu. 

In  this  section  we  presented  a  JTF  ISAR  algorithm  to 
process  data  from  complex  targets  containing  not  only 
scattering  centers,  but  also  other  frequency-dependent 
scattering  mechanisms.  The  adaptive  JTF  ISAR  algo¬ 
rithm  allows  the  enhancement  of  the  ISAR  image  by 
eliminating  nonpoint-scatterer  signals,  thus  leading  to  a 
much  cleaner  ISAR  image.  Moreover,  the  extracted 
nonpoint-scattering  mechanisms  can  be  more  appropri¬ 
ately  displayed  in  an  alternative  feature  space  to  show  tar¬ 
get  resonances  and  frequency  dispersion.  This  is 
accomplished  without  any  loss  in  down-range  resolution. 
In  the  next  section,  the  JTF  processing  will  be  applied  to 
the  cross  range  (or  Doppler  frequency)  instead  of  the 
range  dimension  of  the  radar  image  for  the  purpose  of  re¬ 
focusing  the  image  due  to  complex  target  motion. 

Applying  JTF  Analysis  to  Radar 
Imaging  of  Moving  Targets 

Time-Frequency-Based  Radar  Image  Formation 

To  achieve  high  cross-range  resolution,  SAR  uses  syn¬ 
thetic  array  processing  that  coherently  combines  signals 
obtained  from  sequences  of  small  apertures  at  different 
aspect  angles,  with  respect  to  the  target  to  emulate  the  re¬ 
sult  from  a  large  aperture  [14],  [15].  In  radar  images,  the 
cross-range  resolution  is  obtained  from  the  Doppler  reso¬ 
lution  [13],  [15].  The  differential  Doppler  shifts  of  adja¬ 
cent  scatterers  of  the  target  can  be  observed  in  the  radar 
receiver;  therefore,  the  distribution  of  the  target’s  reflec¬ 
tivity  can  be  measured  by  Doppler  spectra.  The  conven¬ 
tional  method  to  retrieve  Doppler  information  is  the 
Fourier  transform. 

Based  on  the  radar-received  signal  from  a  single  point 
scatterer,  the  received  signal  from  the  entire  target  can  be 
represented  as  the  integration  of  the  returned  signals 
from  all  scatterers  in  the  target  [15] 

S( /)=[[  p(x,y)cxp\-j2Kf—\dxdy 

[  c  }  (9) 

where  r  is  the  range  of  a  point  scatterer  in  the  target.  It  can 
be  expressed  by  the  range  of  the  target’s  center  oi  rotation 
R(f),  the  rotation  angle  of  the  target  0(/“),  and  the  scat¬ 
tered  location  (x,y)  in  the  target  inertial  coordinate  (Fig. 
10).  Because  the  range  R(t)  and  the  rotation  angle  0(r) 


▲  10.  Geometry  of  radar  imaging. 


are  functions  of  time,  the  range  of  a  point  scatterer  in  the 
target  is  also  a  function  of  time 

r  (t)=R(t)+xco$8(t)-ysmQ(t).  (10) 

Then,  the  baseband  signal  in  the  radar  receiver  becomes 

exp i \jjp(x,y)exp{-j2n[xf  (t)]}dxdy 

l  c  (11) 

where  the  frequency  components  are  defined  by 
fx  =— cosG(t) 

7  *  c  (12) 

and 


In  order  to  apply  the  Fourier  transform  properly,  scat¬ 
tered  must  remain  in  their  range  cells  during  the  imaging 
time,  and  their  Doppler  frequency  shifts  must  be  constant 
[17]-[20].  To  obtain  a  focused  Fourier  radar  image, 
range  tracking  must  be  used  to  pull  scattered  back  into 
their  range  cells  and  Doppler  tracking  must  be  used  to 
correct  any  time-varying  phase  change.  Thus,  a  constraint 
has  been  imposed  on  the  Doppler  frequency  spectrum, 
requiring  it  to  be  constant.  Then,  a  clear  radar  image  can 
be  captured  (Fig.  11). 

The  range  tracking  and  Doppler  tracking  form  the  ba¬ 
sis  of  the  standard  procedure  for  focusing  radar  images 
[13],  [15].  If  the  target  is  moving  smoothly,  the  standard 
focusing  is  good  enough  to  generate  a  clear  image  of  the 
target  by  using  the  Fourier  transform.  However,  when  a 
target  exhibits  complex  motion,  such  as  rotation  and  ma¬ 
neuvering,  the  standard  procedure  is  not  sufficient  to 
generate  an  acceptable  image  for  viewing  and  analysis. 
Scattered  can  still  drift  out  of  their  range  cells  and  their 
Doppler  frequency  shifts  can  still  be  time-varying.  Thus, 
the  Doppler  spectrum  obtained  from  the  Fourier  trans¬ 
form  will  be  smeared,  and  the  radar  image  will  be  blurred. 
In  this  case,  more  sophisticated  procedures  for  individual 
scattered,  such  as  polar  reformatting  [16],  are  needed. 
Thus,  each  scatterer  may  remain  in  its  range  cell  and  its 
Doppler  frequency  shift  may  be  constant.  Then,  the  Fou¬ 
rier  transform  can  be  applied  properly  to  reconstruct  a 
clear  image  of  the  target.  However,  to  perform  the  polar 
reformatting,  the  knowledge  required  of  the  initial  kine¬ 
matic  parameters  of  the  target  may  not  be  always  avail¬ 
able.  Some  individual  scattered  may  still  drift  through 


fy  =— sin0(f ). 

(13) 

From  the  radar-received  sig¬ 
nal  S(ft),  if  the  target's  range 
as  a  function  of  time  t  is 
known  exacdy  over  the  imag¬ 
ing  time  duration,  then  the  ex- 
traneous  phase  term 
exp {-j4fR(t)/c}  can  be  re¬ 
moved  by  multiplying  its  con- 
jugate  with  the  received 
_  signal,  i.e.,  G(£,  fy)  =  S(f,t) 
exp{/4/R(t)/r}.  This  is  re¬ 
ferred  to  as  gross  focusing  or 
motion  compensation  [15], 
[18].  Then,  the  target’s  reflec¬ 
tivity  density  function  can  be 
reconstructed  simply  by  tak¬ 
ing  the  inverse  Fourier  trans¬ 
form  of  the  motion 
compensated  signal 

p(x,yi=IFT[G(fx ,/,)]. 

(14) 


their  range  cells  and  their  Doppler  frequencv  shifts  may 


All.  Conventional  radar  imaging  of  moving  targets  with  the  Fourier  transform. 
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A  12.  Time-varying  spectrum  of  radar  data  represented  by  (a) 
the  Fourier  transform  and  (b)  the  time-frequency  transform. 


still  be  time-varying.  Hence  the  resulting  image  can  still 
be  blurred  when  the  Fourier 'transform  is  applied. 

However,  the  restrictions  of  the  Fourier  transform  can 
be  circumvented  if  the  time-frequencv  transform  de¬ 
scribed  in  [1]  is  used  to  replace  the  Fourier  transform 
[48]-[51].  Due  to  the  time- varying  behavior  of  the 
Doppler  frequency  shift,  an  efficient  method  of  solving 
the  problem  of  the  smeared  Fourier  frequencv  spectrum 
and,  hence,  the  blurred  image,  is  to  apply  a 
high- resolution  time-frequencv  transform  to  the  Doppler 
processing.  By  replacing  the  Fourier  transform  with  a 
high-resolution  time-frequency  transform,  the  image 
blurring  caused  by  the  time-varying  Doppler  frequencv 
shifts  can  be  mitigated  without  applying  sophisticated  fo¬ 
cusing  algorithms. 

The  Doppler  frequency  spectrum  of  moving  targets  is 
always  time-varying.  The  relationship  between  the 
time-varying  Doppler  spectrum  and  the  target's  motion 
can  be  described  by  the  translation- induced  Doppler  shift 


and  the  rotation-induced  Doppler  shift  [51] 

2/ 

— [*(-Qsin y  0  -Q.  /rosy  0  )-y(Qcos\\t  0  -D. 2  ts in y  0 )]  ^ 

where  y0  is  the  initial  aspect  angle  of  the  target.  It  is  clear 
that  even  if  the  rotation  rate  Q.  and  the  radial  velocity  of 
the  target  Vr  are  constant,  the  rotation-induced  Doppler 
frequency  shift fdrot  is  still  time- varying.  Other  sources  of 
time-variation  in  the  Doppler  frequency  shift  may  result 
from  uncompensated  phase  errors  due  to  irregularities  in 
the  target's  motion,  the  fluctuation  of  the  tar¬ 
get's  rotation  rate,  fluctuation  in  localizing  the 
rotation  center,  inaccuracy  in  tracking  the 
phase  history  (i.e.,  a  collection  of  phases  of  re¬ 
turned  signals  during  imaging  time),  and  other 
variations  of  the  system  and  the  environment. 
Because  the  residual  phase  errors  may  vary 
with  time,  the  Doppler  frequency  can  also  vary 
with  time. 

The  Fourier  transform  indicates  which  fre¬ 
quency  components  are  contained  in  the  sig¬ 
nal,  but  it  does  not  tell  how  frequencies  change 
with  time.  By  representing  the  time-varying 
Doppler  frequency  spectrum  with  the  Fourier 
transform,  the  Doppler  spectrum  becomes 
smeared  as  m  the  example  shown  in  Fig.  12, 
where  the  Fourier  transform  (a)  and  a 
time-frequency  transform  (b)  are  applied  to  a 
radar  data.  We  can  see  that  the  Fourier  trans¬ 
form  of  the  data  is  actually  the  integral  of  the 
time-frequency  transform  of  the  same  data 
over  its  duration.  This  is  due  to  the  frequency 
marginal  condition  [51]. 

From  the  insight  into  the  image  blurring,  we  can  see 
that  in  order  to  achieve  a  clear  image,  the  time-frequency 
transform  should  be  used  in  place  of  the  Fourier  trans¬ 
form.  The  time-frequency  transform  introduced  in  [1]  is 
an  efficient  way  to  resolve  the  image  blurring  problem 


A  14.  Rotation  motion  represented  by  roll,  pitch  and  yaw. 
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(b)  ISAR  Time-Frequency  Based  Images 


▲  14.  Rotation  motion  represented  by  roll,  pitch ,  and  yaw. 

caused  by  the  time-varying  Doppler’s  behavior  without 
applying  sophisticated  focusing  algorithms  for  individual 
scatterers. 

Fig.  13  illustrates  the  time-frequency-based  radar  im¬ 
age  formation  [50],  [51].  Standard  range  tracking  and 
Doppler  tracking  are  necessary  prior  to  performing  the 
time-frequency  image  formation.  The  difference  between 
the  time-frequency-based  image  formation  and  the  con¬ 
ventional  Fourier-based  image  formation  is  that  the  fo¬ 
cusing  procedure  for  individual  scatterers  is  no  longer 
needed  in  the  time-frequency-based  image  formation  and 
the  Fourier  transform  is  replaced  by  the  time-frequency 
transform  followed  by  time  sampling.  The  Fourier-based 
imaging  approach  generates  only  one  image  frame  from  a 
radar  data  set  G{rm(n )},  where  rm(n)  is  the  m- th  range 
profile  of  Af  range  profiles  and  n  =  However, 

the  time-frequency-based  image  formation  takes 
time-frequency  transforms  at  each  range  cell  and  gener¬ 
ates  an  NxN  Doppler  frequency  (or  Doppler)  and 
dwell-time  (or  time)  distribution.  By  combining  the  M 
Doppler  time  distributions  at  Af-range  cells,  theMXNx 
AT-range-Doppler-time  cube  Q(rmfn,tn)  can  be  formed 


Q(rm  J„Jn  )=TF{G[rm  (»)]} 


(17) 


where  TF  denotes  the  time-frequency  operation  with  re¬ 
spect  to  n.  At  a  particular  time  instant  f-,  only  one  range 
Doppler  image  frame  Q(rm^fn,  tn  =  ts)  can  be  extracted 
from  the  cube.  There  are  a  total  of N  image  frames  avail¬ 
able,  and  every  one  represents  a  full  range-Doppler  image 
at  a  particular  instant.  Therefore,  by  replacing  the  Fourier 
transform  with  the  time-frequency  transform,  a  2-D 
range-Doppler  Fourier  image  frame  becomes  a  3D 
range-DoppIer-time  image  cube.  By  sampling  in  time,  a 
time  sequence  of  2-D  range  Doppler  images  can  be 
viewed  [49]-[51].  Each  individual  time-sampled  frame 
from  the  cube  provides  not  only  a  clear  image  with  supe¬ 
rior  resolution,  but  also  time-varving  properties  from  one 
time  to  another. 


Time-Frequency  Analysis  for 
Moving  Targets 

A  moving  target  always  has  translational  and 
rotational  motions  during  the  imaging  time. 
Assume  the  radar  is  located  at  the  origin  of 
the  (u,v,w)  frame  of  reference,  and  a  target  is 
located  at  a  distance  R(u,v,w,t )  from  the  ra¬ 
dar  at  time  t.  The  inertial  coordinate  of  the 
target  is  the  (x,y,z)  frame.  Thus,  rotational 
motions  of  the  target  can  be  described  by 
roll,  pitch,  and  yaw  [15],  [18].  For  an  air¬ 
craft  heading  along  the  jc-axis,  roll  corre¬ 
sponds  to  a  rotation  0r  about  the.x;-axis,  pitch 
corresponds  to  a  rotation  0p  about  they-axis, 
and  yaw  corresponds  to  a  rotation  0V  about 
the  2-axis  (Fig.  14).  If  the  order  of  rotations 
is  a  roll,  followed  by  a  pitch,  and  finally,  a 
yaw,  then  the  composite  roll,  pitch,  and  yaw  motions  can 
be  represented  by  a  rotation  matrix 

Rot(Qr  ,0 ,  ,0,  )=Rollx  (0r  )Pitcby  (0 ,  )Tawz  (0, ).  (18) 

According  to  this  rotation  matrix,  any  point  scatterer  on 
the  target  will  move  from  its  current  location  to  a  new  lo¬ 
cation.  When  the  target  is  rotating,  the  Doppler  fre¬ 
quency  spectrum  of  the  returned  signal  from  the  target 
becomes  time-varying,  and  the  reconstructed  image  will 
be  blurred.  For  example,  if  a  target  has  pitch  motion  only, 
i.e.,  0r  =  0,  =  0  and  0^,  =  £lp  t ,  where  Clp  is  the  pitch  rate, 
the  rotation  matrix  becomes 


Rot(Qr  ,e,,e,)= 


’  cos0, 

0 

sinG, 

O' 

0 

1 

0 

0 

-sinO, 

0 

COS0, 

0 

0 

0 

0 

1 

(19) 


For  a  point  scatterer  at  x,y,z,  the  Doppler  frequency  shift 
induced  by  the  pitch  motion  becomes 


(20) 


where  we  assume  the  initial  pitch  angle  is  zero.  Even  if  the 
pitch  rate  Qp  is  a  constant,  the  motion-induced  Doppler 
frequency  shift  fd  ^  is  still  time-varying.  Assume  the 
point  scatterer  at  Je  f=  5  m,  2  =  0  and  the  target  rotating 
with  a  pitch  rate  of  £2p  =  0.14  rad/s,  from  (20),  the  Dopp¬ 
ler  bandwidth  is  approximately  6  Hz.  Consequently,  if 
the  Doppler  resolution  is  1  Hz,  corresponding  to  a  imag¬ 
ing  time  of  1  s,  the  motion-induced  Doppler  frequency  of 
the  scatterer  will  be  smeared  over  6  Doppler  resolution 
cells,  and  the  measurement  of  its  true  position  will  be  un¬ 
certain. 

Fig.  15(a)  shows  the  conventional  ISAR  Fourier  im¬ 
age  of  a  target  taken  by  a  X-band  radar  operating  at  9000 
MHz  [52].  Because  the  target  is  maneuvering,  the  Fou¬ 
rier  image  of  the  target  is  still  blurred  even  after  applying 


* i  'iDfu  1  mo 


FfGNAI  PROCFSCINC  MAra 7!NF 


From  the  insight  into  the  image 
blurring,  we  can  see  that  in  order 
to  achieve  a  clear  image,  the 
time-frequency  transform  should 
be  used  in  place  of  the  Fourier 
transform. 

the  standard  procedure  for  focusing.  By  using  the 
time-frequency-based  image  formation,  because  each 
scatterer  has  its  own  range  and  Doppler  frequency  shift  at 
each  time  instant,  without  knowing  any  targets  kine¬ 
matic  parameters  and  resampling  the  data,  a  blurred  Fou¬ 
rier  image  becomes  a  sequence  of  clear  time- frequency 
images.  Fig.  15(b)  shows  one  of  these  time-frequency  im¬ 
ages  in  which  the  nose,  wing-tips,  fuselage,  and  engines 
of  the  aircraft  can  be  seen  very  clearly.  The  blurred  image 
due  to  target  maneuvering  can  be  refocused  without  ap¬ 
plying  sophisticated  autofocusing  or  motion  compensa¬ 
tion  algorithms. 

Summary 

The  Fourier  transform  has  been  widely  used  in  radar  sig¬ 
nal  and  image  processing.  When  the  radar  signals  exhibit 
time-  or  frequency-varying  behavior,  an  analysis  that  can 
represent  the  intensity  or  energy  distribution  of  signals  in 
the  joint  time-frequency  domain  is  most  desirable.  In  this 
article,  we  showed  that  JTF  analysis  is  a  useful  tool  for  im¬ 
proving  radar  signal  and  image  processing  for  time-  and 
frequency-varying  cases.  We  applied  JTF  analysis  to  radar 
backscattering  and  feature  extraction;  we  also  examined 
its  application  to  radar  imaging  of  moving  targets.  Other 
applications  of  time-frequency  analysis  to  radar  can  also 
be  found  in  [53]-[56]. 

Most  methods  of  JTF  analysis  are  non-parametric. 
However,  parametric  or  model-based  methods  of 
time-frequency  analysis,  such  as  adaptive  Gaussian  and 
chirplets,  are  more  suitable  for  radar  signals  and  images 
[57]-[59]. 
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ABSTRACT 

Based  on  the  adaptive  joint  time-frequency  processing  techniques,  a  new  methodology  is  proposed  in  this  paper  to  separate 
the  interference  due  to  fast  rotating  parts  from  the  original  ISAR  image  of  the  target.  The  technique  entails  adaptively 
searching  for  the  linear  chirp  bases  which  best  represent  the  time- frequency  behavior  of  the  signal  and  fully  parameterizing 
the  signal  with  these  basis  functions.  The  signal  components  due  to  the  fast  rotating  part  are  considered  to  be  associated  with 
those  chirp  bases  having  large  displacement  and  slope  parameters,  while  the  signal  components  due  to  the  target  body  motion 
are  represented  by  those  chirp  bases  which  have  relatively  small  displacement  and  slope  parameters.  By  sorting  these  chirp 
bases  according  to  their  slopes  and  displacements,  the  scattering  due  to  the  fast  rotating  part  can  be  separated  from  that  due  to 
the  target  body.  Consequently,  the  image  artifacts  overlapping  with  the  original  image  of  the  target  can  be  well  removed  and 
a  clean  ISAR  image  can  be  produced.  Successful  applications  of  the  algorithm  to  numerically  simulated  and  measurement 
data  show  the  robustness  of  the  algorithm. 

Keywords:  ISAR  imaging,  Doppler  smearing,  fast  rotating  parts,  adaptive  joint  time-frequency  processing 


1.  INTRODUCTION 

It  is  well  known  that  when  rotating  components  exist  on  a  target  such  as  gimbaled  antennas  or  propeller  blades,  image 
artifacts  are  introduced  in  the  Doppler  dimension  of  the  inverse  synthetic  aperture  radar  (ISAR)  image1' 2’ 3.  The  reason  is 
those  parts  are  often  rotating  at  a  speed  much  faster  than  the  target  itself  and  the  Doppler  frequencies  involved  by  them  are 
often  much  higher  than  the  radar  pulse  repetition  rate.  Therefore  serious  aliasing  can  happen  and  result  in  smearing  through 
all  the  cross  range  dimension  of  the  image.  These  smeared  features  oftentimes  overshadow  the  target  geometrical  features 
and  hinder  the  proper  interpretation  of  the  ISAR  image.  To  enhance  such  an  image,  the  most  often  used  approach  is  to  simply 
estimate  the  size  of  the  target  in  the  cross  range  dimension,  keep  the  signal  components  inside  the  target  dimension  and 
discard  all  the  other  unwanted  components  outside  the  confined  region.  By  this  approach  one  can  get  rid  of  most  of  the 
Doppler  smearing  if  the  actual  size  of  the  target  is  known.  However,  in  most  of  the  real  world  ISAR  applications,  the  target  is 
unknown  and  it  is  difficult  to  give  a  close  estimation  of  the  target  size.  An  over  estimation  of  the  target  size  must  be  used  to 
assure  most  of  the  target  features  are  kept  in  the  image,  while  it  may  result  in  poor  elimination  of  the  Doppler  smearing. 
Besides,  even  if  the  size  of  the  target  is  accurately  known,  those  Doppler  smearing  overlapping  with  target  features  can  never 
be  eliminated  by  this  gating  technique,  since  they  are  within  the  geometrical  dimension  of  the  target.  To  overcome  this 
difficulty,  a  technique  is  presented  in  this  paper  to-remove~Such  Doppler  smearing  and  produce  a  clear  ISAR  image  of  the 
target  based  on  the  adaptive  joint  time-frequency  processing. 

The  original  concept  of  adaptive  joint  time-frequency  processing  was  proposed  in  the  signal  processing  society  4  5  and 
applied  to  ISAR  image  processing  in  the  joint  time-frequency  space  for  resonant  scattering  mechanism  extraction  6’  7  and 
target  motion  compensation  8‘  9.  Comparing  to  the  regular  processing  in  either  the  time  or  the  frequency  domain,  more 
insights  into  the  underlying  physical  mechanisms  can  be  gained  in  the  joint  time-frequency  plane  because  it  takes  the 
advantage  of  the  instantaneous  frequency  behaviors  of  the  signal.  The  technique  entails  adaptively  searching  for  the  linear 


chirp  bases  which  best  represent  the  time-frequency  behavior  of  the  signal.  This  is  accomplished  by  projecting  the  signal  onto 
all  possible  chirp  bases  and  finding  the  one  with  the  maximum  projection  value.  After  the  optimal  basis  is  found,  the  signal 
component  associated  with  this  basis  is  subtracted  from  the  original  signal.  By  iterating  this  search  procedure,  the  signal  can 
be  fully  parameterized  with  a  set  of  chirp  basis  functions.  Since  the  Doppler  frequency  due  to  the  rotating  component  is  both 
higher  and  more  rapidly  varying  (in  dwell  time)  than  that  from  the  target  body,  the  signal  components  due  to  the  fast  rotating 
part  are  associated  with  those  chirp  bases  having  large  displacement  and  slope  parameters.  On  the  other  hand,  those  chirp 
bases  which  have  relatively  small  displacement  and  slope  parameters  represent  the  signal  components  due  to  the  target  body 
motion.  By  sorting  these  chirp  bases  according  to  their  slopes  and  displacements,  the  scattering  due  to  the  fast  rotating  part 
can  be  separated  from  that  due  to  the  target  body.  Consequently  a  cleaned  ISAR  image  of  the  target  body  can  be 
reconstructed  by  using  only  those  bases  associated  with  the  body  motion.  Furthermore,  if  the  pulse  repetition  rate  (PRF)  of 
the  radar  is  high  enough,  an  image  of  the  rotating  part  can  also  be  constructed  by  removing  the  body  interference.  We  shall 
present  the  results  of  applying  this  algorithm  to  simulated  data  from  the  radar  scattering  prediction  code  Xpatch. 


2.  ADAPTIVE  JOINT  TIME-FREQUENCY  PROCESSING 

ISAR  imaging  is  a  simple  and  robust  process  for  mapping  the  position  and  magnitude  of  the  target  scattering  features  based 
on  a  point  scatterer  assumption  of  the  target.  With  this  assumption,  the  target  consists  of  a  group  of  point  scatterers. 
Therefore  the  radar  scattering  signal  collected  versus  frequency  and  observation  angle  can  be  described  as: 


E(f,  9)  =  '^0(xi,yi)e-2JkXlCOse-2lky'sine 

1=1 


(1) 


where  0{xx  ,y,j  is  the  amplitude  of  the  ith  scattering  center,  k  is  the  free  space  wave  number.  In  real  radar  measurements  for 
targets  in  flight,  the  looking  angle  is  usually  unknown  but  scaled  according  to  the  observation  dwell  time.  If  the  target  is 
rotating  with  respect  to  the  radar,  the  signal  can  be  written  as: 


i  y 

E(f,tD)  =  ^0(xi,yi)e 


-2  jkXj  cos (OfD  )-2  jky{  sin(0/D ) 


/=! 


(2) 


where  Q  is  the  angular  rotation  velocity  of  the  target.  Note  that  if  the  sinusoid  functions  are  simplified  under  the  small 
observation  angle  approximation,  the  range  and  cross  range  information  can  be  obtained  by  Fourier  transforming  the  signal 
from  the  frequency  and  dwell  time  domain  *.  However,  when  fast  rotating  parts  exist  on  the  target,  eq.(2)  doesn’t  hold  any 
more  and  has  to  be  rewritten  to  include  the  scattering  from  the  parts  that  are  rotating  at  a  different  velocity.  We  limit  the 
discussion  to  the  case  when  only  one  rotation  speed  exists  for  those  parts  other  than  the  target  motion  and  replace  (2)  by  the 
following: 


A 

,)  =  2^0(xi. 


-2jkx,  cos(nbtD)-2jkVj  sin (QbtD) 


i=\ 


^Oix^y^e 


-2jkx{  cos(  Q,ptD)-2  jky{  sin (ClptD ) 


(3) 


where  Nb  is  assumed  to  be  the  number  of  all  the  scatterers  in  the  target  body,  Qband  Qpare  respectively  the  angular  velocities 
of  the  target  body  and  the  fast  rotating  parts.  Usually  £2p  is  much  greater  than  Qb.  In  order  to  image  the  target,  the  observation 
time  is  chosen  to  allow  the  target  to  be  observed  for  a  small  range  of  angles.  While  the  first  term  in  (3)  can  still  be  projected 
to  the  image  of  the  target  via  the  Fourier  transform,  the  second  term  may  result  in  serious  Doppler  smearing  in  cross  range 
domain  and  overshadow  the  target  features  since  the  small  angle  approximation  does  not  hold  any  more.  For  this  reason,  the 
signal  have  to  be  separated  for  different  motion  mechanisms  before  the  regular  Fourier  processing  can  be  applied. 

We  shall  utilize  the  adaptive  joint  time-frequency  processing  technique  to  achieve  this  goal.  The  algorithm  implemented  here 
is  very  similar  to  that  in  our  earlier  work  8.  The  criterion  to  tell  the  two  components  of  the  signal  in  (3)  is  according  to  their 


different  Doppler  frequencies  versus  dwell  time  characteristics.  For  the  component  due  to  target  body  scattering,  the  Doppler 
frequency  is: 

fc  =  2KcQb[ycos(QhtD)  +  xsin(QbtD)]  =  2Kcnb(y  +  xQbtD)  (4) 

which  is  generally  a  linear  function  of  dwell  time  with  very  small  slope  and  displacement  parameters  because  Qb  is  very 
small  and  the  small  angle  approximation  can  be  applied.  On  the  other  hand,  the  Doppler  frequency  due  to  the  fast  rotating 
parts  is  basically  a  sinusoid  function  of  dwell  time  as: 

fo  =2KcClp[ycos(QptD)  +  xsm(QptD)]  (5) 

Note  that  both  the  amplitude  and  frequency  of  this  sinusoid 
are  scaled  by  Qp.  This  means  the  Doppler  frequency  can 
either  be  very  large  or  fast  varying  versus  dwell  time.  Thus 
the  signal  can  also  be  characterized  by  the  linear  chirp  bases 
in  dwell  time-Doppler  frequency  plane  but  with  either  large 
displacement  parameters  or  large  slope  parameters. 

Consequently,  if  the  signal  can  be  fully  parameterized  by 
these  bases,  the  component  due  to  the  fast  rotating  parts  can 
be  easily  separated  from  the  component  due  to  target  body 
motion  by  sorting  them  into  different  categories  according 
to  their  parameters.  This  concept  is  depicted  in  Fig.  1  to  state 
the  different  time-frequency  behaviors  of  the  signal 
components  due  to  the  steady  scattering  centers  on  the 
target  body  and  the  fast  rotating  parts.  To  carry  out  the 
parameterization  and  sorting  process,  we  utilize  the 
adaptive  joint  time-frequency  processing  technique  to  be 
described  as  follows. 

First,  the  chirp  basis  functions  with  linear  characteristics  in 
the  (dwell  time)-(Doppler  frequency)  plane  are  constructed 
with  only  linear  and  quadratic  phase  terms: 

hp(tD)  =  exp[-y2ff(/0ro  +\f,tD2)]  (6) 

where  the  linear  phase  coefficient  f0  is  the  displacement  parameter  in  the  time-frequency  plane  and  the  quadratic  phase 
coefficient  f  is  the  slope  parameter.  We  shall  process  the  Doppler  signal  for  each  range  cell  independently  instead  of  for  each 
frequency,  since  the  Doppler  smearing  due  to  the  rotating  parts  is  well  isolated  in  a  finite  number  of  range  cells.  Thus  we 
search  the  chirp  basis  function  which  best  represents  the  time-frequency  behavior  of  the  signal  for  each  range  cell  by 
maximize  the  projection  value  of  the  signal  to  the  basis: 

\BP\2  =max/o/i  \\Rp(x,tD)h'p(tD)dt\  (7) 

where  R(x,  tD)  is  the  range  profile  which  is  the  1-D  Fourier  transform  of  the  original  signal  E(f,  tD)  from  the  frequency 
domain  to  the  range  domain.  The  search  for  the  first  order  coefficient  can  be  accomplished  by  using  the  FFT  algorithm.  Then 
only  a  one-dimensional  search  is  required  to  find  f.  This  procedure  is  equivalent  to  picking  out  the  strongest  chirp 
component  of  the  signal  in  the  time-frequency  plane  with  the  resolution  of  the  full  Doppler  bandwidth.  The  index  p  denotes 
the  signal  is  in  the  p\h  stage  of  the  iterative  procedure  since  this  process  will  be  repeated  for  the  residue  signal  which  is: 


Doppler  frequency 

Fig.l.  The  time-frequency  characteristics  of  the 
components  due  to  scattering  centers  and 
rotating  parts. 


RP+i  (x,tD)  =  Rp  (x,  tD)-Bp-hp  (tD ) 


(8) 


Therefore ,  this i  searching-extraction  process  will  be  iterated  until  the  energy  of  the  residue  signal  is  smaller  than  a  preset 
threshold.  At  this  stage,  the  signal  has  been  fully  parameterized  as  a  sum  of  all  these  bases.  The  signal  component  due  to  the 
target  body  scattering  can  thus  be  reconstructed  by  summing  all  the  bases  with  small  displacement  and  small  slope 
parameters.  In  this  manner,  a  clean  ISAR  image  of  the  target  body  can  be  obtained  via  the  standard  Fourier  processing. 

3.  NUMERICAL  SIMULATION  AND  EXAMPLES 


Cross  range 


Fig. 2.  The  geometry  of  the  helicopter  model. 


Fig. 3.  The  primitive  ISAR  image  obtained  directly 
via  the  Fourier  transform. 


To  verify  the  algorithm  and  better  explain  the 
advantage  of  the  algorithm,  we  have  simulated  the 
radar  scattering  from  a  helicopter  model  using  an 
electromagnetic  prediction  code  Xpatch  which  is 
based  on  the  shooting  and  bouncing  ray  technique 
and  tested  our  algorithm  using  the  simulation 
data.  The  geometry  of  the  helicopter  model  is 
plotted  in  Fig.2.  There  are  four  blades  on  the  rotor. 
The  length  of  the  blade  is  10  feet.  The  target  body 
is  simply  made  by  connecting  two  conducting 
boxes,  which  are  respectively  6.56x6.56x6.56  feet 
and  9.84x3.28x3.28  feet.  In  the  Xpatch  simulation, 
the  radar  scattering  is  observed  for  128  frequency 
points  from  1.75  GHz  to  2.25  GHz  and  128  points 
in  target  aspect  from  158  degrees  to  170  degrees, 
but  from  0  degree  to  360  degrees  for  the  blade 
aspect.  The  elevation  aspect  angle  of  the  target  is  20 
degrees.  Under  this  parameters,  if  the  rotation 
velocity  of  the  target  body  is  12  rpm,  the  blade  will 
be  spinning  at  a  speed  30  times  of  it,  which  is  360 
rpm.  The  standard  ISAR  image  of  the  target  is 
obtained  via  a  2-D  Fourier  transform  as  shown  in 
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Fig.4.  The  dwell  time-Doppler  frequency  spectrogram  of  the 
signal  at  the  range  cell  60. 


Fig. 3.  The  target  body  scattering  features  are 
confined  to  a  small  extent  around  the  center  of  the 
cross  range,  while  the  scattering  components  due  to 
the  fast  rotating  blades  behave  like  strong  Doppler 
smearing  lines  running  across  the  whole  cross  range 
domain  and  overlapping  with  the  scattering  features 
of  the  target  body.  To  observe  the  (dwell  time)- 
(Doppier  frequency)  behaviors  of  the  signal  with 
Doppler  smearing,  we  choose  the  range  cell  60  and 
use  the  short  time  Fourier  transform  to  display  the 
time  frequency  spectrogram  in  Fig.  4.  Although  the 
resolution  is  bad  due  to  the  short  time  Fourier 
transform,  we  can  still  see  there  are  two  kinds  of 
time-frequency  behaviors  going  on  in  Fig.4.  One  is 
those  components  around  the  zero  Doppler 
frequency  (the  center),  which  behave  like  straight 
lines  with  very  little  slope.  It  represents  the 
scattering  from  the  scattering  centers  on  the  target 
body.  The  other  time-frequency  mechanism  is  those 
lines  going  through  the  entire  Doppler  frequency 
domain.  Instead  of  the  sinusoid  like  time-frequency 
behavior  mentioned  in  the  last  section,  these  lines 
are  the  results  of  serious  aliasing  effect  in  the 
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Fig.5.  The  ISAR  image  obtained  by  gating  out  the 
Doppler  component  outside  the  estimated 
target  region. 


Doppler  frequency  Domain.  The  aliasing  occurs  whenever  the  highest  Doppler  frequency  induced  by  the  fast  rotating  parts  is 
greater  than  the  radar  pulse  repetition  rate  (PRF),  which  causes  serious  aliasing  in  Doppler  frequency  domain.  Even  though 
the  sinusoid  behavior  can  no  longer  be  seen,  the  linear  chirp  description  of  the  time-frequency  behavior  due  to  the  rotating 
parts  is  still  effective.  The  chirp  basis  functions,  however,  have  to  be  defined  with  very  large  slope  parameters  to  well 
represent  the  aliasing  effect. 


As  far  as  the  ISAR  image  in  Fig.3  is  concerned,  the 
standard  way  to  process  it  is  to  estimate  the  target 
dimension  and  the  corresponding  Doppler  frequency 
extent,  and  then  gate  out  the  unwanted  Doppler 
frequency  components.  For  this  example,  we  choose 
a  Doppler  frequency  window  which  allows  about 
twice  of  the  target  dimension  and  zoom-in  on  this 
region.  The  resulting  image  is  plotted  as  Fig.5.  Note 
that  little  is  done  to  eliminate  the  Doppler  smearing 
that  overlaps  with  the  target  features.  Now  the 
adaptive  joint  time-frequency  processing  is  applied 
to  the  original  signal.  First  we  Fourier  transform  the 
original  frequency-aspect  data  to  the  range-aspect 
domain.  Then  for  each  range  cell,  the  signal  is 
projected  to  all  possible  bases  until  the  optimal 
chirp  basis  is  found  and  extracted.  The  procedure  is 
repeated  until  the  residue  energy  is  less  than  2%  of 
the  energy  of  the  strongest  range  cell.  Usually,  it 
takes  less  than  200  terms  to  parameterize  the  signal 
with  the  residue  energy  less  than  such  a  threshold.  It 
has  been  noted  that  in  the  search  the  slope 
parameters  of  the  bases  should  be  large  enough  to 
represent  the  dramatic  phase  variation  of  the 
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Fig.6.  The  ISAR  image  after  applying  the  adaptive 
joint  time- frequency  processing  to  remove 
the  Doppler  smearing. 


Doppler  signals  due  to  the  fast  rotating  parts.  After  the  parameterization  is  done,  the  basis  functions  are  divided  into  two 
categories  according  to  their  slope  and  displacement  parameters.  Here  we  choose  an  appropriate  small  number  as  the 
threshold  of  slope.  The  threshold  of  displacement  is  chosen  the  same  as  in  the  standard  approach.  Only  those  bases  with  both 
the  slope  and  the  displacement  lower  than  the  thresholds  are  kept  as  the  components  due  to  the  scattering  from  the  target 


body-  Summing  up  these  bases  and  Fourier  transforming  it  to  the  image  domain,  a  clean  ISAR  image  is  obtained  in  Fig.6.  We 
observe  that  most  of  the  Doppler  smearing  is  eliminated.  For  a  better  assessment  of  the  Doppler  smearing  removal,  the  ISAR 
image  simulated  for  the  target  without  the  rotating  blades  is  also  generated  in  Fig.7.  From  the  comparison,  it  is  seen  that  the 
algorithm  has  successfully  kept  almost  all  the  scattering  features  due  to  the  target  body,  but  has  removed  most  of  the  Doppler 
smearing  interference  due  to  the  rotating  parts  from  the  image. 


4.  IMAGE  RECONSTRUCTION  FOR  THE  ROTATING  PARTS 


In  the  above  adaptive  joint  time-frequency  processing,  the  scattering  component  due  to  the  rotating  parts  can  also  be  isolated 
in  the  process.  It  consists  of  the  bases  with  either  the 
slope  parameters  or  the  displacement  parameters 
larger  than  the  thresholds.  If  the  signal  is  seriously 
aliased,  it  is  difficult  to  gather  any  useful 
information  on  the  rotating  parts.  Nevertheless,  if 
the  PRF  of  the  radar  is  high  enough  to  avoid 
aliasing,  we  should  be  able  to  generate  an  image  for 
the  rotating  parts  as  well.  However,  just  as  in  the 
case  of  the  Doppler  signal  from  the  rotating  parts 
causing  artifacts  in  the  image  of  the  target  body,  the  „ 

scattering  of  the  target  body  can  also  cause  artifacts  c 

in  the  image  of  the  blades.  For  example,  we  simulate  a! 
the  radar  scattering  of  the  same  target  again  using 
Xpatch  for  the  same  frequency  sampling  but  128 
observation  angles  from  0  degree  to  12  degrees  for 
the  blades.  In  this  range,  the  body  rotates  only  0.4 
degrees.  Although  the  blades  can  be  imaged  clearly, 
the  poor  resolution  for  the  body  compresses  the 
body  image  into  zero  cross  range  as  shown  in  Fig.8. 

■■F  Cro.  Si;.nge  (dB) 


In  conjunction  with  the  information  provided  by  the 
parameterization,  it  is  also  possible  to  suppress  this 
body  interference.  Note  that  the  scattering 
component  due  to  the  target  body  is  parameterized 
as  a  sum  of  chirp  basis  functions  versus  dwell  time. 
The  parameterization  model  expression  is  still 
effective  even  for  the  shorter  observation  period. 
The  scattered  fields  due  to  the  target  body  can  thus 
be  interpolated  to  the  dense  sampling  used  to 
observe  the  blades  and  be  subtracted  from  the  total 
scattering  fields.  The  ISAR  image  of  the  blades  with 
the  body  interference  suppressed  is  obtained  as 
Fig-9.  Comparing  the  results  with  the  ISAR  image 
simulated  with  the  blades  only  in  Fig.  10,  the  body 
line  interference  is  much  weaker  but  not  completely 
eliminated.  The  reason  could  be  the  effects  due  to 
the  multiple  scattering  mechanism  simulated  by 
Xpatch.  in  which  the  total  scattered  fields  can  no 
longer  be  written  as  the  sum  of  the  scattered  fields 
due  to  the  target  body  and  due  to  the  blades. 


Fig.7.  The  ISAR  image  of  the  target  body  simulated 
without  the  rotating  blades. 
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Fig.8.  The  ISAR  image  of  the  blades  with  the  body 
interference. 


5.  SUMMARY 


In  this  paper  we  have  applied  adaptive  joint  time- 
frequency  ideas  to  separate  two  different  scattering 
mechanisms  respectively  due  to  the  target  body  and 
the  fast  rotating  parts  on  the  target.  This  is 
accomplished  by  taking  the  advantage  of  their 
different  (dwell  time)-(Doppler  frequency) 
behaviors.  After  adaptively  searching  for  the  linear 
chirp  bases  which  best  represent  the  time-frequency 
behavior  of  the  signal  and  fully  parameterizing  the 
signal  with  these  basis  functions,  we  can  sort  these 
bases  into  the  signal  component  due  to  the  fast 
rotating  parts  and  the  component  due  to  target  body 
motion  according  to  their  slope  and  displacement 
parameters.  The  Doppler  smearing  interference 
caused  by  the  fast  rotating  parts  can  be  removed 
from  the  original  ISAR  image  by  keeping  only  those 
bases  associated  with  target  body  scattering.  The 
algorithm  has  been  tested  for  numerical  simulated 
radar  scattering  data  for  a  helicopter  model.  The 
ISAR  image  obtained  by  removing  the  Doppler 
smearing  from  the  original  image  agrees  well  with 
the  image  simulated  for  the  target  without  the 
rotating  blades  on.  Furthermore,  the  algorithm  has 
also  been  extended  to  image  the  fast  rotating  blades 
with  body  interference  removed  for  radar  with  high 
PRF.  This  algorithm  has  been  applied  to  the 
measurernent  data  for  real  targets.  These  results  will 
be  presented  in  the  talk. 
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Fig.9.  The  ISAR  image  of  the  blades  with  the  body 
interference  removed. 
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Fig.  10.  The  ISAR  image  simulated  with  the  blades  only. 
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ABSTRACT 

The  application  of  joint  time-frequency  techniques  for  the  analysis  of  electromagnetic  backscattered  data  is  reviewed.  In  the 
joint  time-frequency  feature  space,  discrete  time  events  such  as  scattering  centers,  discrete  frequency  events  such  as  target 
resonances,  and  dispersive  mechanisms  due  to  surface  waves  and  guided  modes  can  be  simultaneously  displayed.  We 
discuss  the  various  joint  time-frequency  representations  including  the  short-time  Fourier  transform,  wavelet  transform, 
Wigner-Ville  distribution,  windowed  super-resolution  algorithms  and  the  adaptive  spectrogram.  Emphasis  is  placed  on  how 
these  algorithms  can  be  used  to  represent  with  good  resolution  the  scattering  phenomenology  in  electromagnetic  data.  We 
highlight  an  application  of  joint  time-frequency  processing  for  radar  image  enhancement  and  feature  extraction.  It  is  shown 
that  by  applying  joint  time-frequency  processing  to  the  conventional  inverse  synthetic  aperture  radar  (ISA  R)  imagery,  it  is 
possible  to  remove  non-point  scattering  features  in  the  image,  leading  to  a  cleaned  image  containing  only  physically 
meaningful  point  scatterers.  The  extracted  frequency-dependent  mechanisms  can  be  displayed  in  an  alternative  feature  space  to 
facilitate  target  identification. 

Keywords:  joint  time- frequency  analysis,  electromagnetic  scattering,  feature  extraction,  radar  imaging 


1.  INTRODUCTION 

The  usefulness  of  the  joint  time-frequency  analysis  of  signals  has  long  been  recognized  in  the  signal  processing  arenas.1,2 
Recently,  joint  time-frequency  methods  have  also  been  applied  to  electromagnetic  data  with  good  success.  The  joint  time- 
frequency  representation  of  a  signal  is  a  two-dimensional  phase  space  that  facilitates  the  visualization  and  interpretation  of 
complex  electromagnetic  wave  phenomenology.  In  this  feature  space,  discrete  time  events  such  as  scattering  centers,  discrete 
frequency  events  such  as  resonances,  and  dispersive  mechanisms  due  to  surface  waves  and  guided  modes  can  be 
simultaneously  displayed.  This  can  oftentimes  lead  to  more  insights  into  the  complex  electromagnetic  wave  propagation  and 
scattering  mechanisms  than  what  is  available  in  the  traditional  time  or  frequency  domain  alone.  The  objectives  of  this  paper 
are  to  provide  an  overview  of  the  various  joint  time-frequency  techniques  and  to  show  how  they  can  be  applied  to  extract 
meaningful  phenomenology-based  features  from  electromagnetic  data.  We  will  discuss  various  time- frequency  representations 
of  signals  including  the  short-time  Fourier  transform,  wavelet  transform,  Wigner-Ville  distribution,  windowed  super¬ 
resolution  algorithms  and  the  adaptive  spectrogram.  Emphasis  will  be  placed  on  how  electromagnetic  phenomenology  is 
manifested  in  the  joint  time-frequency  plane  and  how  these  features  can  be  interpreted  and  extracted.  We  will  highlight  an 
application  of  joint  time-frequency  processing  for  radar  image  enhancement  and  feature  extraction. 

The  standard  tool  in  generating  the  joint  time-frequency  representation  of  a  signal  is  the  short-time  Fourier  transform 
(STFT).  First  introduced  by  Gabor  in  1946,  it  is  basically  a  sliding  window  Fourier  transform  in  time.  First,  let  us  define 
the  Founer  transform  of  a  signal  f(t): 

F(co)  =  J  f(t)  e'J®1  dt  (1) 


By  taking  the  Fourier  transform  of  the  windowed  time  data  as  the  window  is  moved  in  time,  a  two-dimensional  time- 
frequency  image,  or  the  spectrogram,  is  generated: 


STFT  (t,  £1)  =  J  f(t)  g(t-x)  e'J^  dt  (2) 

Equation  (2)  is  very  similar  to  the  Fourier  transform  defined  in  (1),  except  for  the  presence  of  a  frequency  window  function 
g(t).  The  definition  of  the  STFT  can  also  be  expressed  in  the  time  domain  by  manipulating  (3a),  with  the  result: 

STFT  (x,  Q)  =  e'jflT  J  F(o>)  G(f2-co)  e*™  do)  (3) 

Here  G(oi)  is  the  inverse  Fourier  transform  of  g(t).  The  dual  relationship  between  (2)  and  (3)  is  apparent,  i.e.,  the  time- 
frequency  representation  can  be  obtained  through  either  a  moving  window  in  time  or  a  moving  window  in  frequency.  In 
addition,  we  make  the  following  observations:  (i)  Time  signals  with  duration  shorter  than  the  duration  of  the  window  will 
tend  to  get  smeared  out,  i.e.,  the  resolution  in  the  time  domain  is  limited  by  the  width  of  the  window  g(t).  Similarly  the 
resolution  in  the  frequency  domain  is  limited  by  the  width  of  the  frequency  window  G(co).  (ii)  The  window  width  in  time 
and  the  window  width  in  frequencyare  inversely  proportional  to  each  other  by  the  uncertainty  principle.  Therefore,  good 
resolution  in  time  (small  time  window)  necessarily  implies  poor  resolution  in  frequency  (large  frequency  window),  (iii)  The 
window  width  in  each  domain  remains  fixed  as  it  is  translated.  This  results  in  a  fixed  resolution  across  the  entire  time- 
frequency  plane. 
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Fig.  1 .  Electromagnetic  mechanisms  are  manifested  in  the  joint  time-frequency  image  as  distinct  features. 

The  spectrogram  provides  information  on  the  frequency  content  of  the  signal  at  different  time  instances.  In 
electromagnetic  applications,  the  joint  time-frequency  image  can  be  used  to  simultaneously  display  the  different  wave 
phenomenology  due  to  discrete  events  of  time,  discrete  events  of  frequency  and  dispersive  mechanisms.'5'4  Shown  in  Fig.  1 
are  the  time-frequency  features  of  some  common  mechanisms.  A  discrete  event  of  time,  for  example,  could  be  due  to  wave 
scattering  from  a  spatially  localized  "scattering  center"  on  a  structure.  It  shows  up  as  a  vertical  line  (Fig.  1(a))  in  the  image 
since  it  occurs  at  a  particular  time  instance  but  over  all  frequencies.  A  resonance,  which  is  a  discrete  event  of  frequency, 
shows  up  as  a  horizontal  line  (Fig.  1(b)).  Dispersive  phenomena,  on  the  other  hand,  are  characterized  by  slanted  curves  in 
the  time-frequency  image.  For  example,  surface  wave  mechanisms  due  to  material  coatings  are  characterized  by  curves  with  a 
positive  slope  in  the  time-frequency  image  (Fig.  1(c)).  As  the  frequency  is  increased,  the  surfacewave  becomes  more  tightly 
bound  to  the  material  layer,  and  consequently,  propagates  at  a  slower  velocity  and  results  in  a  longer  time  delay.  Another 
type  of  dispersion  is  due  to  waveguide  structures.  These  "structural  dispersion"  mechanisms  are  characterized  by  curves  with 


a  negative  slope  in  the  time- frequency  image  (Fig.  1(d)).  This  can  be  explained  by  the  fact  that  as  the  frequency  is  lowered, 
one  approaches  the  cutoff  frequency  of  a  waveguide  mode  and  the  propagation  delay  associated  with  this  mechanism 
approaches  infinity.  As  we  shall  see  in  this  paper,  all  of  the  above  mentioned  phenomena  have  been  observed  in  a  wide 
variety  of  structures,  from  simulation  data  on  canonical  structures  to  measurement  data  on  complex  platforms.  For  radar 
applications,  for  instance,  the  time-frequency  representation  is  particularly  effectivefor  identifying  scattering  mechanisms  in 
targets  containing  "sub-skinline"  structures  such  as  inlet  ducts,  antenna  windows  and  material  coatings. 


2.  JOINT  TIME-FREQUENCY  ALGORITHMS  AND  SCATTERING  PHENOMENOLOGY 


The  insights  gained  in  the  time-frequency  plane  come  at  the  price  of  resolution.  As  we  have  already  discussed,  the  time- 
frequency  image  generated  by  the  STFT  is  limited  in  resolution  by  the  extent  of  the  sliding  window  function.  Smaller  time 
window  results  in  better  time  resolution,  but  leads  to  worse  frequency  resolution,  and  vice  versa.  To  overcome  the 
resolution  limit  of  the  STFT,  a  wealth  of  alternative  time- frequency  representations  have  been  developed  in  the  signal 
processing  community.  Below,  we  shall  introduce  four  classes  of  algorithms  which  have  been  utilized  to  analyze 
electromagnetic  data,  namely,  continuous  wavelet  transform,  Wigner-Ville  distribution,  windowed  superresolution 
algorithms  and  adaptive  techniques.  While  this  list  is  by  no  means  complete,  these  four  classes  of  algorithms  provide  a 
flavor  of  the  different  approaches  to  the  resolution  issue. 

2.1  Continuous  Wavelet  Transform.  Contrary  to  the  fixed  resolution  of  the  STFT,  the  wavelet  transform  can  be  defined 
to  achieve  variable  resolution  in  one  domain  (either  time  or  frequency)and  multi-resolution  in  the  other  domain.  We  will 
now  introduce  the  continuous  wavelet  transform  (C  WT),  or  the  so-called  scalogram,  of  a  frequency  signal  F(co): 


CWT(t,fl)=  —  f  F((0)  t1/2  d(t) 

2kj 


(4) 


H(co)  is  usually  termed  the  "mother  wavelet"  in  wavelet  theory.  Essentially,  (4)  can  be  interpreted  as  a  decomposition  of  the 
frequency  signal  F(co)  into  a  family  of  shifted  and  dilated  wavelets  H(x(co-h)).  The  wavelet  basis  function  H(x(co-£2))  has 
variable  w  *dth  according  to  x  at  each  frequency  £2.  H(x(co-£2))  is  wide  for  small  x  and  narrow  for  large  x.  By  shif  ir*^  H(co) 
with  a  fixed  scale  parameter  t,  the  1/x-scale  mechanisms  in  the  frequency  response  F(co)  can  be  extracted  and  localized. 
Alternatively,  by  dilating  H(co)  at  a  fixed  £2,  all  of  the  multi-scale  events  of  F(co)  at  £2  can  be  analyzed  according  to  the  scale 
parameter  x.  This  is  the  so-called  "multi-resolution"  property  of  the  wavelet  transform  and  is  an  advantage  over  the  STFT 
for  analyzing  multi-scale  signals.  Like  the  STFT,  the  wavelet  transform  can  also  be  carried  out  on  the  inverse  Fourier 
transform  f(t)  of  the  original  frequency  signal  F(co): 

CWT  =  J  f(t)  t-1/2  h  (-t  /  x)  e‘jQt  dt  (5) 


where  h(t)  is  the  Fourier  transform  of  H(co).  Since  (4)  is  essentially  the  Fourier  transform  of  [f(t)x~l/2h(-t/x)],  it  is  the 
preferred  numerical  implementation  of  the  wavelet  transform  through  the  use  of  the  fast  Fourier  transform  algorithm  for  each 
value  of  x.  By  comparing  Equations  (2)  and  (5),  we  note  that  h(t)  is  similar  to  the  window  function  in  the  STFT.  However, 
h(t)  must  satisfy  the  "admissibility  condition"  in  wavelet  theory,  viz.,  h(t=0)  =  0.  To  satisfy  this  condition,  h(t)  is  usually 
chosen  to  be  a  shifted  window  function  with  its  center  at  to-  By  changing  x,  this  window  function  is  shifted  by  xto  and  the 
width  of  the  window  is  dilated  by  the  scale  factor  x.  The  ratio  between  the  window  width  and  the  window  center  (or  the  Q- 
factor  of  the  window  function)  remains  fixed  for  all  times.  This  is  in  contrast  to  the  STFT  where  the  window  width  does  not 
change  as  it  is  being  shifted. 

It  is  worthwhile  to  point  out  here  that  the  definition  of  the  wavelet  transform  presented  above  in  its  time  and 
frequency  forms  is  exactly  the  complement  of  its  usual  definition  used  in  time-series  signal  analysis.  For  electromagnetic 
applications,  it  has  been  found  that  the  property  of  the  wavelet  transform  we  are  usually  interested  in  is  its  multi-resolution 
capability  in  the  frequency  domain  and  its  variable  resolution  capability*  in  the  time  domain.  The  multi-resolution  capability 
is  ideally  suited  for  analyzing  frequency-domain  electromagnetic  data  which  consist  of  both  discrete  time  events  (of  large 
extent  in  frequency)  and  discrete  frequency  events  (of  small  extent  in  frequency). 

As  an  example,  the  time-frequency  representation  of  the  backscattering  data  from  an  open-ended  waveguide  duct  is 
considered.5  The  duct  is  an  open-ended  circular  waveguide  with  a  diameter  of  1.75  in.  A  flat  conducting  termination  exists 
2  ft.  inside  the  waveguide  (Fig.  2(a)).  To  generate  the  backscattering  data,  the  radar  cross  section  of  this  target  is  first 
computed  in  the  frequency  domain.  The  time-domain  response  is  then  obtained  by  Fourier  transforming  the  band-limited 


frequency  data  (from  2  to  18  GHz).  Fig.  2(b)  shows  the  spectrogram  of  the  backscattering  data  for  the  H-polarization  at  45° 
off-normal  incidence  using  the  STFT.  In  performing  the  STFT,  a  2-GHz  Kaiser-Bessel  window  in  the  frequency  domain  is 
used.  Also  plotted  along  the  two  axes  are  the  time-domain  and  the  frequency-domainresponses.  It  is  apparent  that  the 
scattering  features  are  much  better  resolved  in  the  time-frequency  domain  than  in  either  the  time  or  the  frequency  domain 
alone.  Both  the  non-dispersive  rim  diffraction  and  the  mode  spectra  due  to  multiple  propagating  modes  in  the  circular 
waveguide  can  be  clearly  identified.  The  mode  spectra  are  in  fact  dispersion  curves  of  the  waveguide  modes  since  the  phase 
velocity  of  each  mode  is  proportional  to  the  travel  time. 
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Fig.  2.  Joint  time-frequency  images  of  the  backscattered  data  from  an  open-ended  waveguide  duct 
obtained  via  STFT  and  continuous  wavelet  transform. 


Due  to  the  fixed  resolution  of  the  STFT,  the  scattering  features  in  Fig.  2(a)  are  smeared  out  in  the  time-frequency 
plane.  This  problem  is  overcome  by  using  the  wavelet  transform  which  provides  a  much  better  representation  of  the 
scattering  features  in  the  time-frequency  plane,  as  shown  in  Fig.  2(c).  The  wavelet  transform  is  implemented  using  (5)  with 
the  aid  of  the  FFT.  The  function  h(t)  is  chosen  to  be  a  two-sided  Kaiser-Bessel  window  with  a  Q-factorof  0.3.  The  t=0 
reference  ofh(t)  is  located  midway  between  the  time  events  from  the  rim  diffraction  and  interior  contribution  (at  x=2.05  ns). 
The  variable  time  resolution  of  the  wavelet  transform  allows  sharper  time  resolution  to  be  achieved  during  the  early-time 
response  and  sharper  frequency  resolution  (coarser  time  resolution)  to  be  achieved  during  the  late-time  response.  Thus, 
wavelet  transform  provides  good  resolution  in  identifying  the  scattering  centers  and  resolving  the  resonant  phenomena  of  the 
target  while  adequately  describing  the  dispersive  scattering  mechanisms  in  the  intermediate-time  region. 

2.2  Wigner-Ville  Distribution  and  Time-Frequency  Distribution  Series.  While  the  STFT  and  the  wavelet  transform  are 
based  on  linear  transformations,  another  class  of  time-frequency  distributions  can  be  obtained  from  the  quadratic,  or  power 
spectrum,  point  of  view.  The  most  basic  of  these  is  called  the  Wigner-Ville  distribution  (WVD),  first  developed  in  quantum 
mechanics.  The  WVD  of  a  signal  f(t)  is  defined  as: 

WVD(t,Q)=  J  f(x  +  ^  )  )  e'Jat  dt  (6) 


The  WVD  distribution  is  always  real  and  gives  information  on  how  the  power  spectrum  of  the  signal  changes  as  a  function  of 
time.  Although  the  WVD  gives  nearly  the  best  resolution  among  all  the  time-frequency  techniques,  its  main  drawback 
comes  from  cross-term  interferenceproblem.  Simpiv  put,  the  WVD  of  the  sum  of  two  signals  is  not  the  sum  of  their 


WVD‘s.  Therefore,  if  a  signal  contains  more  than  one  component  in  the  joint  time-frequency  plane,  its  WVD  will  contain 
cross  terms  which  occur  halfway  between  each  pair  of  the  auto  terms.  The  magnitude  of  these  oscillatory  cross  terms  can  be 
twice  as  large  as  the  auto  terms  and  yet  they  do  not  have  any  physical  meaning.  This  drawback  severely  limits  the 
usefulness  of  the  WVD  in  its  original  form.  A  number  of  techniques  can  be  used  to  alleviate  this  problem.  We  shall 
introduce  one  recently  proposed  by  Qian  and  Chen.6  They  suggested  that  if  the  WVD  can  be  decomposed  into  a  sum  of 
localized  and  symmetric  functions,  it  may  be  possible  to  suppress  cross-term  interference  by  selecting  only  the  low  order 
harmonics.  This  is  accomplished  by  first  decomposing  the  original  signal  into  the  Gabor  expansion: 


W  ^  ^  ^m,n  ^m,n  (0 

m  n 

where 

Vn  (0  =  (*a2)-0-25  exp  +  JnAt  }  (8) 

are  time-shifted  and  frequency-modulated  Gaussian  basis  functions.  In  the  above  expression,  m  and  T  are  respectively  the 
time  sampling  index  and  time  sampling  interval,  while  n  and  A  are  the  sampling  index  and  sampling  interval  in  frequency. 
By  taking  the  WVD  of  both  sides  of  (7),  it  becomes: 

WVD  (x,Q)  =  £  X  Cm,nC*m',n’  WVD^.  (t,Q)  (9) 

mn  m’  n' 

where  WVDh  h.  denotes  the  WVD  between  any  pair  of  basis  functions  and  is  available  in  closed  form.  Next,  the  above 
expression  can  be  regrouped  based  on  the  “interaction  distance” 


D  =  |  m  -  m'  |  +  |  n  -  n'  |  (10) 

between  the  pairs  of  bases  (m,  n)  and  (m\  n').  This  results  in  what  is  termed  the  "time-frequency  distribution  series"  (TFDS, 
also  called  the  Gabor  spectrogram): 

WVD(x.Q)  =  X  I  Cm,n  I2  WVDy,.  (*.Q)  (D=0  terms)  (11) 

mn 


+  X  X  Cm,n  cV„.  WVDhh.(x,Q)  (D=l  terms) 

mn  m'n1 


+  X  X  cm,n  cVn'  WVDh  h.  (x,£2)  (D=2  terms) 

mn  m’n’ 


Clearly,  if  we  take  all  the  terms  in  the  series  (D=°°),  the  right  hand  side  of  (1 1)  converges  to  the  WVD  of  the  original  signal. 
This  yields  the  best  resolution  but  is  plagued  by  cross-term  interference.  At  the  other  extreme,  if  we  take  only  the  self¬ 
interaction  terms  in  the  series  (D=0),  the  resulting  right-hand  side  is  equivalent  to  the  STFT  of  the  signal  using  a  Gaussian 
window  function.  It  has  no  cross-term  interference  problem  but  has  the  worst  resolution.  Therefore,  as  the  order  D  increases, 
we  gain  in  resolution  but  pay  a  price  in  cross-term  interference.  It  is  often  possible  to  balance  the  resolution  against  cross¬ 
term  interference  by  adjusting  the  order  D.  The  optimal  value  for  D  was  reported  to  be  around  2  to  4. 

As  an  example,  we  consider  the  scattering  from  a  dielectric  coated  wire.  The  coated  wire  is  a  12"  section  ot  a 
coaxial  cable  (RG-41/U  with  a=0.037"  and  b= 0.095")  with  its  outer  conductor  removed.  The  dielectric  material  is  Teflon 


with  permittivity  er= 2.1.  The  scattering  data  are  collected  by  both  numerical  simulation  and  by  chamber  measurement.  The 
spectrograms  for  the  case  of  60°  incidence  from  broadside  is  shown  in  Fig.  3(a).  we  observe  a  slight  tilting  of  the  vertical 
lines  in  the  spectrograms,  especially  in  the  late-time  returns.  This  signifies  the  presence  of  dispersive  phenomena  which  are 
due  to  the  surface  wave  mechanism  in  the  dielectric  coating  known  as  the  Goubau  mode.  At  frequencies  well  above  DC,  this 
mode  is  tightly  bound  to  the  dielectric  and  the  group  velocity  of  the  wave  approaches  the  slow  dielectric  velocity.  As  the 
frequency  approaches  DC,  the  wave  velocity  approaches  that  of  free  space  and  exhibits  a  shorter  propagation  delay.  Therefore, 
in  the  time-frequency  plane,  the  Goubau  mode  phenomena  show  up  as  slanted  curves  with  a  positive  slope.  The  TFDS  and 
the  WVD  results  are  also  given  in  Figs.  3(b)  and  3(c),  respectively.  As  we  can  see  from  Figs.  3(c),  the  WVD  of  the  data  is 
so  contaminated  by  cross-term  interference  it  is  essentially  useless.  The  TFDS  results  shown  in  Fig.  3(b)  is  a  good 
compromise  between  resolution  and  cross-term  interference.  The  returns  that  can  not  be  distinguished  from  each  other  in  the 
spectrogram  are  now  easily  identified  in  the  joint  time-frequency  plane  by  using  the  TFDS  of  order  2. 
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Fig.  3.  Joint  time-frequency  representations  of  the  backscattered  data  from  a  dielectric-coated  wire  obtained  via  STFT,  time- 
frequency  distribution  series  with  N=2,  and  Wigner-Ville  distribution. 

2.3  Windowed  Superresolution  Algorithms.  In  either  the  STFT  or  the  wavelet  transform,  the  resolution  in  the  time- 
frequency  display  is  limited  by  the  extent  of  the  sliding  window  function  as  opposed  to  the  full  bandwidth  of  the  signal.  The 
processing  of  the  data  within  each  time  or  frequency  window  using  a  super-resolution  technique  such  as  ESPRIT,  MUSIC  or 
matrix-pencil  algorithm  therefore  appears  quite  attractive.8,9  Super-resolved  parameterization  retains  the  advantage  of 
simultaneous  time-frequency  display  while  overcoming  the  resolution  issue.  However,  additional  processing  is  needed  to 
fully  parameterize  the  data,  especially  when  dispersive  mechanisms  are  present.  Furthermore,  the  robustness  of  the 
algorithms  to  noise  needs  to  be  carefully  considered.  We  will  describe  a  simple  windowed  superresolution  procedure  based 
on  Prony's  method  for  achieving  parameter  estimation  of  both  scattering  centers  and  natural  resonances  in  the  time-frequency 
plane. 


In  the  windowed  time-frequency  super-resolution  procedure,  Prony's  extraction  is  first  applied  in  the  frequency 
domain  to  locate  discrete  time  events.  Prony's  method  will  fit  the  raw  data  to  the  following  model: 

-jcotm 

F(0))=2.Ame  (12) 

m=l 


where  the  tm's  are  the  unknown  locations  of  the  discrete  time  events,  the  A,n’s  are  the  unknown  complex  strengths  of  these 
events,  and  M  is  the  number  of  them  to  be  found.  It  is  clear  that  applying  Prony's  method  to  the  entire  frequency  data  will 


yield  a  poor  fit  if  the  raw  data  contains  discrete  frequency  events  such  as  resonances.  To  circumvent  this  problem,  we 
repeated  apply  Prony’s  method  to  many  small  windows  of  the  raw  data.  Prony’s  method  will  be  successful  for  most  of  the 
window  locations  and  will  fail  only  when  a  window  coincides  with  a  resonant  peak.  By  repeatedly  sliding  the  window  along 
in  frequency  and  re-applying  Prony’s  method  to  the  raw  data,  we  are  able  to  identify  as  true  locations  those  values  of  tm 
which  most  frequently  occur.  A  weighted  least  squares  fit  of  the  values  of  Am  for  each  discrete  time  event  is  used  to  construct 
a  smooth  functional  form  of  Am.  An  important  benefit  of  this  global  functional  form,  Am(co),  is  that  it  allows  us  to  go  back 
and  interpolate  Am  at  those  frequencies  where  Prony’s  method  originally  failed. 

Provided  that  the  discrete  time  events  have  been  properly  located  and  accounted  for,  the  remaining  data  will  consist 
solely  of  a  series  of  resonant  peaks  after  the  frequency  domain  extraction.  To  extract  the  natural  resonance  information,  this 
remaining  data  is  inverse  Fourier  transformed  to  the  time  domain.  The  sliding  window  Prony's  procedure  is  then  applied  in 
the  time  domain  to  fit  the  complex- valued  data  to  a  model  which  is  the  dual  of  (12): 

VI  -joint 

f(t)=XBne  (13) 

n=  1 


in  which  the  con’s  are  the  unknown  resonance  frequencies,  and  the  Bn’s  are  their  corresponding  strengths.  By  tracking  the 
behavior  of  each  Bn  with  respect  to  time,  other  parameters  associated  with  the  resonance,  such  as  attenuation  factor  (an)  and 
turn-on  time  (Tq)  are  extracted: 

N  +jCOn(t-Tn)  -an(t-Tn) 

f(t)  =  2.  bn  e  e  U(t— Xn)  (14) 

n=l 

in  which  the  bn’s  are  strengths  of  the  resonances  at  turn-on  and  u^-Tn)  is  the  unit  step  function. 
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Fig.  4.  Time-frequency  images  of  the  backscattered  data  from  a  conducting  strip  with  an  open  cavity 
obtained  via  STFT  and  windowed  superresoiuhon  algorithm. 


As  an  example,  let  us  consider  a  perfectly  conducting  strip  containing  a  small  open  cavity  located  in  the  middle  of 
the  strip  shown  in  Fig.  4(a).  Backscattering  data  from  50  MHz  to  18  GHz  were  simulated  using  a  method  of  moments 
solution  for  the  E-poiarized  incidence  at  25*  from  edge-on.  Shown  in  Fig.  4(b)  is  the  spectrogram  of  the  scattered  data 
obtained  via  the  STFT  with  a  1-GHz  frequency  window.  The  three  vertical  lines  correspond  to  the  scattering  centers  from 
the  left  edge,  the  cavity  exterior  and  the  right  fin,  and  the  five  visible  horizontal  lines  correspond  to  the  cavity  resonances.  It 
is  possible  to  extract  qualitative  information  such  as  the  approximate  frequency  behaviors  of  the  individual  scattering  centers 
and  the  resonance  Q’s  from  the  STFT  image.  However,  due  to  the  large  number  of  features  contained  in  the  backscattered 
data,  the  image  is  blurry,  making  it  difficult  to  resolve  fine  details.  Fig.  4(c)  shows  a  time-frequency  plot  of  the 
parameterizedbackscattereddata.  Because  we  have  completely  parameterized  the  data  via  a  super-resolution  technique,  the 
sharpness  of  the  image  is  not  constrained  by  the  well-known  Fourier  limit  as  is  the  case  for  the  STFT  image.  Our  image  can 
be  of  nearly  infinite  sharpness,  so  we  have  chosen  for  each  mechanism  to  appear  as  either  a  horizontal  or  vertical  line  exactly 
one  pixel  in  width.  The  intensities  of  the  three  vertical  lines  show  that  the  three  scattering  centers  are  of  differing  strengths 
and  have  differentfrequency  behaviors.  The  high-Q  resonances  can  be  seen  to  be  of  much  longer  duration  than  the  low-Q 
resonances. 

2.4  Adaptive  Gaussian  Representation.  While  the  use  of  wavelet  is  a  step  towards  variable  resolution  in  the  time- 
frequency  plane,  it  is  rather  rigid  in  its  particular  form  of  the  time- frequency  grid.  Flexible  resolution  in  the  time-frequency 
plane  to  accommodate  components  of  the  signal  with  differentresolution  is  therefore  highly  desirable.  Recently,  a  signal 
representation  scheme  that  uses  adaptive  normalized  Gaussian  functions  was  introduced/0  (A  similar  algorithm  called 
matching  pursuit  was  developed  independently  at  around  the  same  time  by  Mallat  and  Zhang11).  Unlike  the  STFT  and 
wavelet  decomposition,  the  time  and  frequency  resolution  as  well  as  the  time-frequency  centers  are  adjusted  to  best  match  the 
signal.  The  objective  of  this  method  is  to  expand  a  signal  f(t)  in  terms  of  normalized  Gaussian  functions  hp(t)  with  an 
adjustable  standard  deviation  cp  and  a  time-frequency  center  (tp,  fp): 


oo 

f(t)  =  Z  Bp  hp  (t) 

P-l 

where 

2  -0.25  f  (t  -  U)2  l  /  x 
hp  (t)  =  ( 7tcTp  )  exp[ - J  exp(j  2jcfpt) 

2CTP 


(15) 


(16) 


Note  that  the  modulated  Gaussian  basis  has  a  dual  form  in  its  Fourier  transform  representation: 

Hp  (0  =  (1/2jictp)2]  0  25  exp[ - -  '  p)  2  ]  exp(-j  2n(f-  fp)  tp) 

2(  1  /2tco"p) 


(17) 


Therefore,  these  basis  functions  have  a  time-frequency  extent  given  by,  respectively,  ap  and  (l/2rcap).  The  coefficients  Bp  are 
found  one  at  a  time  by  an  iterative  procedure.  One  begins  at  the  stage  p=l  and  chooses  the  parameters  <jp,  tp  and  fp  such  that 
hp(t)  is  most  “similar”  to  f(t),  that  is: 


max 


(>tp,fp 


h  (t)  dt 


(18) 


where  fo(t)  =f(t).  For  p  >1,  fp(t)  is  the  remainder  after  the  orthogonal  projection  of  fp_|(t)  onto  hp(t)  has  been  removed  from 
the  signal: 


fpW  =  fp-l  W  -  Bp(t)hp(t) 


(19) 


This  procedure  is  iterated  to  generate  as  many  coefficients  as  needed  to  accurately  represent  the  original  signal. 


Several  comments  can  be  made  about  the  adaptive  Gaussian  representation.  First,  it  can  be  shown  that  the  norm  of 
the  residue  monotonically  decreases  and  converges  to  zero,  therefore  adding  a  new  term  in  the  series  does  not  affectthe 
previously  selected  parameters.  Second,  because  this  representation  is  adaptive,  it  will  generally  be  concentrated  in  a  very 
small  subspace.  As  a  result,  we  can  use  a  finite  summation  of  the  terms  in  (15)  to  approximate  the  signal,  with  a  residual 
error  as  small  as  one  wishes.  Also,  since  random  noise  in  general  is  distributed  uniformly  in  the  entire  space,  this  subspace 
representation  actually  increases  the  signal-to-noise  ratio.  Finally,  the  major  difficulty  in  implementing  this  algorithm  is  the 
determination  of  the  optimal  elementary  function  at  each  stage.  One  implementation  strategy  is  to  start  with  a  large  ap  and 
scan  the  data  in  frequency  and  time  fora  peak.  Then  we  divide  op  by  two  and  find  the  new  peak.  We  continue  this 
procedure  until  the  standard  deviation  is  small  enough  and  then  select  the  highest  peak  and  extract  the  residual  using  (19).  It 
should  be  pointed  out  that  the  fast  Fourier  transform  is  used  during  the  search  procedure  to  obtain  the  coefficientsfor  all  the 
frequency  centers  at  once,  speeding  up  the  search  that  would  otherwise  be  very  time  consuming. 

The  result  of  applying  the  adaptive  Gaussian  extraction  can  be  effectively  displayed  in  the  time-frequency  plane  using 
the  so-called  adaptive  spectrogram  (ADS): 

ft  t 

ADS  (t,  f)  =  2  X  |Bpf2  exp[-  — §—  -  (2jc)2  o2p  (f- fp)  2  ]  (20) 

P  aP 

This  representation  was  obtained  by  calculating  the  Wigner-Ville  distribution  of  (15)  and  then  deleting  the  cross  terms.  It  can 
be  shown  that  the  energy  contained  in  the  ADS  is  identical  to  the  energy  contained  in  the  signal.  So  it  can  be  considered  as  a 
signal  energy  distribution  in  the  time- frequency  domain.  It  is  also  non-negative,  cross-term  interference  free,  and  of  high 
resolution. 
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Fig.  5.  Adaptive  spectrogram  of  the  backscattered  data  from  a  conducting  strip 
with  an  open  cavity  (see  Fig.  4(a)). 

To  show  an  example  of  the  adaptive  spectrogram,  we  used  the  same  structure  used  in  Fig.  4(a).  We  applied  the 
algorithm  described  above  and  obtained  an  approximation  of  the  signal  using  50  terms  in  expression  (1)  with  a  residual 
energy  of  0.1%.  The  ADS  obtained  is  shown  in  Fig.  5.  We  can  clearly  see  and  locate  with  high  resolution  the  three 
scattering  centers  and  the  resonances.  We  notice  that  the  high-Q  resonances  appear  as  very  thin  horizontal  lines,  while  the 
low-Q  resonances  appear  as  thicker  lines  as  expected.  We  can  even  see  that  the  second  and  the  third  scattering  center 
corresponding  respectively  to  the  cavity  exterior  and  the  right  fin  are  in  reality  multiple  scattering  centers. 


3.  AN  APPLICATION  -  RADAR  IMAGE  ENHANCEMENT  AND  FEATURE  EXTRACTION 

We  will  now  describe  an  applications  of  joint  time-frequency  processing  to  electromagnetic  backscattered  data  in  the  area  of 
radar  imaaing.12  Inverse  synthetic  aperture  radar  (ISAR)  imaging  has  long  been  used  by  the  microwave  radar  community  for 
object  diagnostic  and  target  identification  applications.  ISAR  is  a  simple  and  very  robust  process  for  mapping  the  position 


and  magnitude  of  the  point  scatterers  on  an  object  from  multi-frequency,  multi-aspect  scattered  field  data.  However,  complex 
targets  often  contain  other  scattering  phenomena  such  as  resonances  and  dispersive  mechanisms,  which  do  not  behave  like 
point  scatterers.  One  important  example  is  the  scattering  from  the  engine  inlet/exhaust  duct  on  aircraft.  It  is  a  dominant 
contributor  to  the  overall  scattering  from  the  target,  yet  its  waveguide-like  structure  and  the  associated  frequency-dependent 
scattering  mechanisms  make  it  a  non-point  scattering  feature.  When  processed  and  displayed  by  the  conventional  ISAR 
algorithm,  the  inlet  return  results  in  an  image  feature  which  (i)  is  not  well-focused,  (ii)  is  not  related  to  the  spatial  location  of 
the  scatterer,  and  (iii)  can  often  obscure  other  important  point  features  on  the  target.  Therefore,  it  would  be  useful  to  develop 
an  algorithm  to  automatically  remove  these  features  from  the  ISAR  image,  thus  leading  to  a  cleaned  ISAR  image  containing 
only  physically  meaningful  point  scatterers.  The  extracted  inlet  features,  when  displayed  in  a  more  meaningful  feature  space, 
can  be  used  to  identify  target  resonances  and  cutoff  phenomena. 

In  this  section,  we  present  a  new  processing  technique,  the  joint  time-frequency  ISAR,  that  combines  the  idea  of 
standard  ISAR  with  the  joint  time-frequency  technique  to  accomplish  the  above  objective.  The  conceptual  idea  behind  the 
joint  time-frequency  ISAR  algorithm  is  to  apply  joint  time-frequency  processing  to  the  range  (or  time)  axis  of  the 
conventional  (range)-(cross  range)  ISAR  image  to  gain  an  additional  frequency  dimension.  The  result  is  a  three-dimensional 
(range)-(crossrange)-(frequency)matrix,witheach(range)-(cross  range)  slice  of  this  matrix  representing  an  ISAR  image  at  a 
particular  frequency.  Consequently,  by  examining  how  the  ISAR  image  varies  with  frequency,  we  can  distinguish  the 
frequency-independentscatteringmechanisms  from  the  ffequency-dependentones.  In  the  actual  implementation  of  the  joint 
time-frequency  ISAR,  the  choice  of  the  joint  time-frequency  processing  engine  is  critical.  Our  choice  is  the  adaptive 
spectrogram  described  in  Section  2.4. 


Fig.  6.  (a)  The  VFY-2 18  model,  (b)  Its  standard  ISAR  image  obtained  for  f  ==  8-16  GHz  and  a  40-degree  angular  window 
centered  at  30  degrees  from  nose-on.  (c)  Enhanced  ISAR  image  and  frequency-aspect  display  obtained  by  applying  the 
Adaptive  Gaussian  Representation  to  the  ISAR  image  in  Fig.  6(b).  The  inlet  cloud  has  been  removed  in  the  original  ISAR 
image,  (d)  The  extracted  resonant  features  of  the  inlet  are  shown  in  the  frequency-aspect  plane. 


The  adaptive  spectrogram  has  two  distinct  advantages  over  conventional  time-frequency  techniques  such  as  the 
short-time  Fourier  transform.  First,  it  is  a  parameterization  procedure  that  results  in  very  high  time-frequency  resolution. 
More  importantly  for  our  application,  the  adaptive  spectrogram  allows  us  to  automatically  distinguish  the  frequency- 
dependent  events  from  the  frequency- independent  ones  through  the  extent  of  the  basis  functions.  From  expression  (16)  it  can 
be  seen  that  scattering  centers,  i.e.,  signals  with  very  narrow  length  in  time,  will  be  well  represented  by  basis  functions  with 
very  small  crp.  Frequency  resonances,  on  the  other  hand,  will  be  better  depicted  by  large  crp.  Therefore,  if  we  reconstruct  the 

1SAR  image  using  only  those  Gaussian  bases  with  small  variances,  a  much  cleaner  image  can  be  obtained  showing  only  the 
scattering  centers.  The  remaining  mechanisms,  i.e.,  those  related  to  the  large  variance  Gaussians,  are  more  meaningful  to 
view  in  a  dual  frequency-aspect  display,  where  resonances  and  other  frequency-dependent  mechanisms  can  be  better  identified. 

The  algorithm  is  demonstrated  using  the  chamber  measurement  data  of  a  1:30  scale  model  Lockheed  VFY-218 
airplane  provided  by  the  EMCC  (Electromagnetic  Code  Consortium).12  The  airplane  has  two  long  engine  inlet  ducts,  as 
show  in  Fig.  6(a),  which  are  rectangular  at  the  open  ends  but  merge  together  into  one  circular  section  before  reaching  a  single¬ 
compressor  face.  As  we  can  clearly  see  in  the  conventional  ISAR  image  of  Fig.  6(b)  for  the  vertical  polarization  at  20°  near 
nose-on,  the  large  cloud  outside  of  the  airframe  structure  is  the  inlet  return.  Fig.  6(c)  shows  the  enhanced  ISAR  image  of 
Fig.  6(b),  obtained  by  applying  the  above  joint  time-frequency  ISAR  algorithm  and  keeping  only  the  small  variance 
Gaussians.  We  see  that  only  the  scattering  center  part  of  the  original  signal  remains  in  the  image,  as  expected.  Notice  that  the 
strong  return  due  to  engine  inlet  has  been  removed,  but  the  scattering  from  the  tail  fin  remains.  Fig.  6(d)  shows  the 
frequency-aspect  display  of  the  high  variance  Gaussians.  A  number  of  equispaced  vertical  lines  can  be  seen  between  10.5  and 
13.5  GHz.  Given  the  dimension  of  the  rectangular  inlet  opening,  we  estimate  that  these  frequencies  correspond 
approximately  to  the  second  cutoff  frequency  of  the  waveguide-like  inlet.  We  have  also  run  our  algorithm  on  the  simulated 
scattering  data,  generated  using  a  numerical  Maxwell’s  solver,  for  a  plate-waveguide  configuration  to  verify  this  claim.  The 
frequency  resonances  displayed  indeed  correspond  very  closely  to  the  theoretical  cutoff  frequencies  of  the  waveguide. 

In  this  section  we  presented  a  joint  time-frequency  ISAR  algorithm  to  process  data  from  complex  targets  containing 
not  only  scattering  centers  but  also  other  frequency-dependentscattering  mechanisms.  The  adaptive  joint  time-frequency 
ISAR  algorithm  allows  the  enhancement  of  the  ISAR  image  by  eliminating  non-point  scatterer  signals,  thus  leading  to  a 
much  cleaner  ISAR  image.  It  also  provides  information  on  the  extracted  frequency-dependent  mechanisms  such  as  resonances 
and  frequency  dispersions.  This  is  accomplished  without  any  loss  in  resolution. 
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ABSTRACT 

A  new  methodology  based  on  adaptive  joint  time-frequency  processing  is  proposed  to  separate  the  interference  due  to  fast 
rotating  parts  from  the  original  ISAR  image  of  the  target.  The  technique  entails  adaptively  searching  for  the  linear  chirp  bases 
which  best  represent  the  time-frequency  behavior  of  the  signal  and  fully  parameterizing  the  signal  with  these  basis  functions. 
The  signal  components  due  to  the  fast  rotating  part  are  considered  to  be  associated  with  those  chirp  bases  having  large 
displacement  and  slope  parameters,  while  the  signal  components  due  to  the  target  body  motion  are  represented  by  those  chirp 
bases  having  relatively  small  displacement  and  slope  parameters.  By  sorting  these  chirp  bases  according  to  their  slopes  and 
displacements,  the  scattering  due  to  the  fast  rotating  part  can  be  separated  from  that  due  to  the  target  body.  Consequently,  the 
image  artifacts  overlapping  with  the  original  image  of  the  target  can  be  removed  and  a  clean  ISAR  image  can  be  produced. 
Furthermore,  useful  rotation  rate  information  contained  in  the  Doppler  signal  can  be  extracted.  Successful  applications  of  the 
algorithm  to  numerically  simulated  and  measurement  data  show  the  robustness  of  the  algorithm. 

Keywords:  ISAR  imaging,  Doppler  smearing,  fast  rotating  parts,  adaptive  joint  time-frequency  processing 


1.  INTRODUCTION 

It  is  well  known  that  when  rotating  components  exist  on  a  target  such  as  gimbaled  antennas  or  propeller  blades,  image 
artifacts  are  introduced  in  the  Doppler  dimension  of  the  inverse  synthetic  aperture  radar  (ISAR)  image1, 2* 3.  The  reason  is  that 
those  parts  are  often  rotating  at  a  speed  much  faster  than  the  target  itself  and  the  Doppler  frequencies  generated  by  them  are 
much  higher  than  the  radar  pulse  repetition  rate.  Therefore  serious  aliasing  can  result  in  smearing  through  all  the  cross  range 
dimension  of  the  image.  These  smeared  features  oftentimes  overshadow  the  target  geometrical  features  and  hinder  the  proper 
interpretation  of  the  ISAR  image.  To  enhance  such  an  image,  the  most  often  used  approach  is  to  simply  estimate  the  size  of 
the  target  in  the  cross  range  dimension,  keep  the  signal  components  inside  the  target  dimension  and  discard  all  the  other 
unwanted  components  outside  the  confined  region.  This  approach  can  remove  most  of  the  Doppler  smearing  if  the  actual  size 
of  the  target  is  known.  However,  in  most  of  the  real  world  ISAR  applications,  the  target  is  unknown  and  it  is  difficult  to  give 
a  close  estimation  of  the  target  size.  An  over  estimation  of  the  target  size  must  be  used  to  assure  that  most  of  the  target 
features  are  kept  in  the  image,  thus  resulting  in  poor  elimination  of  the  Doppler  smearing.  Furthermore,  even  if  the  size  of  the 
target  is  accurately  known,  those  Doppler  smearing  overlapping  with  the  target  features  can  never  be  eliminated  by  this 
gating  technique,  since  they  are  within  the  geometrical  dimension  of  the  target. 

Instead  of  treating  the  Doppler  signal  induced  by  the  rotating  parts  as  noise  in  the  image,  it  can  also  be  considered  as  a  type 
of  radar  signature,  since  it  includes  some  unique  information  of  the  target.  For  example,  jet  engine  modulation  (JEM)  of  radar 
return  has  already  been  used  for  target  identification  purposes  4.  In  reference  4,  the  autocorrelation  sequence  of  the  time 
domain  signal  has  been  used  for  period  detection  of  engine  blades.  However,  we  found  that  when  the  target  body  return  is 
strong  and  overlaps  with  the  rotating  part  return  in  the  Doppler  frequency  spectrum,  the  period  detection  algorithm  is  no 
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extraction. 


The  original  concept  of  adaptive  joint  time-frequency  processing  was  proposed 

been  applied  to  ISAR  image Vising  techniques  in  either  the  time  of  the  frequency 
and  target  motion  compensation  •  Compared  to  m  v  s  time-frequency  plane  since  the 

domain,  more  insights  into  the  underlying  physical  mechanism  g  .  ^  t  application  entails  adaptively 

instantaneous  frequency  behaviors  of  the  s.gnal  of  the  signal.  This  is  accomplished  by 

searching  for  the  linear  chirp  bases  which  best  rePr“eJl  1  one  with  toe  maximum  projection  value.  After  the  optimal 
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procedure,  the  signal  can  be  fully  parametenzed  with i  a  set  ofehup  ***  ftom  "  m„t  bodyi  the  s,g„al 
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components  due  to  the  fast  rotating  par.  are  assumed  to  slope  parameters  are 
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2.  ADAPTIVE  JOINT  TIME-FREQUENCY  PROCESSING 

ISAR  imaging  is  a  “d  rob«t  p^for  mppmg  the >  pompon 

frequency  anS  observation  angle  can  be  described  as: 


E(f,d)  =  ^0(Xi,yi)e 


-2  jkXj  cos  8-2  jtyj  sin  8 


(i) 


at  a  constant  rate  with  respect  to  the  radar,  the  signal  can  be  written  as. 


E(f,tD)  =  ^0(xi,yi)e 


-2  jkx,  cos(Qf0 )  -  2  y'fcv,  sin(£2;0) 


(2) 


twre  n  is  the  angular  rotation  velocity  of  the  targe,.  Note  that  if  the  sinusoid  functions  ^ 

observation  angle  approximation,  the  range  and  cross  range  in  ormation  ca”  e  h  (2)  no  longer  hold  any 

from  the  frequency  and  dwell  time  domain  '.  However,  when  fast  rotating  parts  exist  on  th t^get  eq.(* Men  g  ^  ^  ^ 
more  and  has  to  be  rewritten  to  include  the  scattering  from  the  parts  that  are  rotating  at  a  differe 
dtassl^me  case  when  only  on.  rotation  speed  exists  fo,  those  parts  other  than  the  targe,  motton  and  replace  (2)  by 

following: 


(3) 


i=i 


E(f,tD)  =  XO(x,,}',)^2Mcos<at'o)‘2Ms'n<n‘'o)  + 

^  %  -2  jkXjCOsiClptpyljkyjSiniClptQ) 

i=Nh+\ 


where  Nbis  assumed  to  be  the  number  of  all  the  scatterers  in  the  target  body,  Qband  £2pare  respectively  the  angular  velocities 
of  the  target  body  and  the  fast  rotating  parts.  Usually  Opis  much  greater  than  £2b.  In  order  to  image  the  target,  the  observation 
time  is  chosen  to  allow  the  target  to  be  observed  for  a  small  range  of  angles.  While  the  first  term  m  (3)  can  still  be  proj 
to  the  image  of  the  target  via  the  Fourier  transform,  the  second  term  results  in  serious  Doppler  smearing  in  the  cross  '•ange 
domain  and  can  overshadow  the  target  features.  For  this  reason,  the  signal  should  be  separated  for  different  motion 
mechanisms  before  the  regular  Fourier  processing  can  be  applied. 

We  shall  utilize  the  adaptive  joint  time-frequency  processing  technique  to  achieve  this  goal.  The  dpritfim  implemented  here 
is  very  similar  to  that  in  our  earlier  work 8.  The  criterion  to  distinguish  the  two  components  of  the  signal  m  (3)  is  a cord mg  to 
their  different  Doppler  frequencies  versus  dwell  time  characteristics.  For  the  component  due  to  target  body  scattering,  th 

Doppler  frequency  is: 

fZ=2KeQb[ycos(&btD)  +  xsm(ClbtD)]  =  2KlQb(y  +  xQbtD)  (4) 

It  is  generally  a  linear  function  of  dwell  time  with  very  small  slope  and  displacement  parameters  since  Qb  is  ^  smtdland 
the  small  angle  approximation  can  be  applied.  On  the  other  hand,  the  Doppler  frequency  due  to  the  fast  rotating  parts 
basically  a  sinusoid  function  of  dwell  time: 


fl  =  2KcQp[ycos(QptD)  +  xsin(ClptD)] 


(5) 


Note  that  both  the  amplitude  and  frequency  of  this  sinusoid 
are  scaled  by  Op.  This  means  the  Doppler  frequency  can 
either  be  very  large  or  fast  varying  versus  dwell  time.  Thus 
the  signal  can  also  be  characterized  by  the  linear  chirp  bases 
in  dwell  time-Doppler  frequency  plane  but  with  either  large 
displacement  parameters  or  large  slope  parameters. 
Consequently,  if  the  signal  can  be  fully  parameterized  by 
these  bases,  the  component  due  to  the  fast  rotating  parts  can 
be  easily  separated  from  the  component  due  to  target  body 
motion  by  sorting  them  into  different  categories  according 
to  their  parameters.  This  concept  is  depicted  in  Fig.  1  which 
displays  the  different  time-frequency  behaviors  of  the  signal 
components  due  to  the  steady  scattering  centers  on  the 
target  body  and  the  fast  rotating  parts.  To  carry  out  the 
parameterization  and  sorting  process,  we  utilize  the 
adaptive  joint  time-frequency  processing  technique  to  be 
described  next. 

First,  the  chirp  basis  functions  with  linear  characteristics  in 
the  (dwell  time)-(Doppler  frequency)  plane  are  constructed 
with  linear  and  quadratic  phase  terms: 


Doppler  frequency 

Fig.l.  The  time-frequency  characteristics  of  the 
components  due  to  scattering  centers  and 
rotating  parts. 


h  (tD)  =  zxv[-j2n(f0tD  +i/,f/>-)] 


(6) 
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search  the  chirp  basis  function  which  best  represents  the  time-frequency  behavior  of  the  signal  for  each  range  cell  by 
maximize  the  projection  value  of  the  signal  onto  the  basis: 

\Bp\2  =max/o/i  \jRp(x,tD)h‘p(tD)dt\  (7) 

where  R(x,  to)  is  the  range  profile  set  obtained  by  1-D  Fourier  transforming  the  original  signal  E(f,  to)  from  the  frequency 
domain  to  the  range  domain.  The  search  for  the  first-order  coefficient  can  be  accomplished  by  using  the  FFT  algorithm.  Then 
only  a  one-dimensional  search  is  required  to  find  f,.  This  procedure  is  equivalent  to  picking  out  the  strongest  chirp 
component  of  the  signal  in  the  time-frequency  plane  with  the  resolution  of  the  full  Doppler  bandwidth.  The  index  p  denotes 
that  the  signal  is  at  the  pth  stage  of  the  iterative  procedure.  This  process  is  repeated  for  the  residue  signal  which  is: 

Rp«(xJD)  =  Rp(x,tD)-Bp-hp(tD)  (8) 

This  search-extract  process  will  be  iterated  until  the  energy  of  the  residue  signal  is  smaller  than  a  preset  threshold.  At  this 
stage,  the  signal  has  been  fully  parameterized  as  a  sum  of  all  these  bases.  The  signal  component  due  to  the  target  body 
scattering  can  thus  be  reconstructed  by  summing  all  the  bases  with  small  displacement  and  small  slope  parameters.  In  this 
manner,  a  clean  ISAR  image  of  the  target  body  can  be  obtained  via  the  standard  Fourier  processing. 


3.  NUMERICAL  SIMULATION  AND  EXAMPLES 


To  test  the  algorithm,  we  have  simulated  the  radar  scattering  from  a  helicopter  model  using  an  electromagnetic  prediction 
code  Xpatch  which  is  based  on  the  shooting  and  bouncing  ray  technique  n.  The  geometry  of  the  helicopter  model  is  shown  in 
Fig.2.  There  are  four  blades  on  the  rotor.  The  length  of  the  blade  is  10  feet.  The  target  body  is  simply  made  by  connecting 
two  conducting  boxes,  which  are  respectively  6.56x6.56x6.56  feet  and  9.84x3.28x3.28  feet.  In  the  Xpatch  simulation,  the 
radar  scattering  is  observed  for  128  frequency  points  from  1.75  GHz  to  2.25  GHz,  and  128  points  in  target  aspect  from  158 
degrees  to  170  degree;.  'Hie  elevation  angle  of  the  target  is  20  degrees.  Furthermore,  the  blade  position  is  articulated  between 
0  and  360  degrees  during  the  128  angular  looks.  The  parameters  used  assume  that  if  the  rotation  velocity  of  the  target  body  is 
12  rpm,  the  blade  will  be  spinning  at  a  speed  of  30  times  that,  which  is  360  rpm.  The  standard  ISAR  image  of  the  target 


Rg.2.  Geometry  of  the  helicopter  CAD  model. 


Fig.3.  Original  ISAR  image  obtained  directly  via 
the  Fourier  transform. 
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Dwell  Time 


Fig.4.  The  (dwell  time)  -  (Doppler  frequency)  spectrogram  of 
the  signal  at  the  range  cell  53. 

is  obtained  via  a  2-D  Fourier  transform  as  shown  in  Fig.3. 
The  target  body  scattering  features  are  confined  to  a  small 
extent  around  the  center  of  the  cross  range,  while  the 
scattering  components  due  to  the  fast  rotating  blades 
exhibit  strong  Doppler  smearing  lines  running  across  the 
whole  cross  range  domain  and  overlapping  with  the 
scattering  features  of  the  target  body.  To  observe  the  (dwell 
time)-(Doppler  frequency)  behaviors  of  the  signal  with 
Doppler  smearing,  we  choose  the  range  cell  60  and  use  the 
short  time  Fourier  transform  to  display  the  time-frequency 
spectrogram  in  Fig.  4.  Although  the  resolution  is  poor  due 
to  the  short-time  Fourier  transform,  we  can  see  there  are 
two  kinds  of  time-frequency  behaviors  in  Fig.4.  One 
component  is  located  around  the  zero  Doppler  frequency 
(the  center),  which  behaves  like  vertical  straight  lines  with 
very  little  slope.  It  represents  the  scattering  from  the 
scattering  centers  on  the  target  body.  The  other  time- 
frequency  mechanism  is  those  lines  going  through  the 
entire  Doppler  frequency  domain.  Instead  of  the  sinusoid 
like  time-frequency  behavior  mentioned  in  the  last  section, 
these  lines  are  the  results  of  serious  aliasing  effect  in  the 
Doppler  frequency  domain.  The  aliasing  occurs  whenever 
the  highest  Doppler  frequency  induced  by  the  fast  rotating 


Fig.5.  The  ISAR  image  obtained  by  gating  out  the 
Doppler  component  outside  the  estimated 
target  region. 


Cross  range 


Fig.6.  The  ISAR  image  after  applying  the  adaptive 
joint  time-frequency  processing  to  remove 
the  Doppler  smearing. 
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have  to  be  defined  with  very  large  slope  parameters  to  well  represent  the  aliasing  effect. 
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transform  the  original  frequency-aspect  data  to  the 
range-aspect  domain.  Then  for  each  range  cell,  the 
signal  is  projected  to  all  possible  bases  until  the  optimal 
chirp  basis  is  found  and  extracted.  The  procedure  is 
repeated  until  the  residue  energy  is  less  than  2%  of  the 
energy  of  the  strongest  range  cell.  Usually,  it  takes  less 
than  200  terms  to  parameterize  the  signal  with  the 
residue  energy  less  than  such  a  threshold.  It  has  been 
noted  that  in  the  search  the  slope  parameters  of  the  bases 
should  be  large  enough  to  represent  the  dramatic  phase 
variation  of  the  Doppler  signals  due  to  the  fast  rotating 
parts.  After  the  parameterization  is  done,  the  basis 
functions  are  divided  into  two  categories  according  to 
their  slope  and  displacement  parameters.  Here  we 
choose  an  appropriate  small  number  as  the  threshold  of 
slope.  The  threshold  of  displacement  is  chosen  the  same 
as  in  the  standard  approach.  Only  those  bases  with  both 
the  slope  and  the  displacement  lower  than  the  thresholds 
are  kept  as  the  components  due  to  the  scattering  from 
the  target  body.  Summing  up  these  bases  and  Fourier 
transforming  it  to  the  image  domain,  we  arrive  at  a  clean 
ISAR  image  Fig.6.  We  observe  that  most  of  the  Doppler 

smearing  is  eliminated.  To  better  assess  the  effectiveness  of  removing  the  Doppler  smearing,  the  ISAR  image  is  also 
simulated  for  the  target  without  the  rotating  blades  and  is  shown  in  Fig.7.  From  the  comparison,  it  is  seen  that  the  algorithm 
has  successfully  kept  almost  all  the  scattering  features  due  to  the  target  body,  but  has  removed  most  of  the  Doppler  smearing 
interference  due  to  the  rotating  parts  from  the  image.  We  have  also  tested  the  algorithm  on  measured  ISAR  data  with  good 
success. 


(dB) 


Fig.7.  The  ISAR  image  of  the  target  body  simulated 
without  the  rotating  blades. 


4.  DOPPLER  INFORMATION  EXTRACTION 

In  the  above  AJTF  processing,  the  scattering  component  due  to  the  rotating  parts  can  also  be  isolated  in  the  process.  It 
consists  of  the  bases  with  either  the  slope  parameters  or  the  displacement  parameters  larger  than  the  thresholds.  Because  the 
radar  PRF  is  too  low  to  capture  the  Doppler  frequency  induced  by  the  fast  rotation,  the  signal  is  seriously  aliased.  However, 
the  periodicity  of  the  signal  envelope  due  to  the  periodic  motion  of  the  blades  is  still  available.  Usually  the  observation  time 
period  is  long  enough  to  observe  several  periods  of  the  signal  envelope.  Therefore,  a  period  detection  algorithm  can  be 
developed  to  extract  the  rotation  rate  of  the  rotating  parts. 

The  period  tp  can  be  detected  by  observing  the  amplitude  peaks  in  the  autocorrelation  function  versus  dwell  time,  which  is  : 

f(x,T)=\  j  Rp(x,t0)»  Rp(x,tD-r)dtD  I  (9) 

The  peaks  will  occur  at  X  =  rt  •  t  ,  provided  the  signal  is  purely  from  the  rotating  parts,  where  n  =0,  1,2, . For  multiple 

blade  rotor,  tp  is  equal  to  the  ratio  of  the  rotation  period  and  the  number  of  blades.  Note  the  prerequisite  of  the  algorithm  is 
the  periodicity  of  the  signal  envelope,  which  means  the  signal  should  be  due  only  to  the  rotating  parts  to  achieve  a  robust 
detection.  Any  signal  interference  from  the  target  body  will  degrade  the  performance  of  the  algorithm.  Thus  AJTF  processing 
is  used  here  as  an  effective  tool  to  remove  the  body  interference.. To  show  the  advantage  of  the  adaptive  joint  time-frequency 
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processing,  we  applied  the  algorithm  to  a  set  of  real 
measurement  result.  The  measurement  is  from  a  helicopter.  The 
autocorrelation  function  of  the  original  signal  versus  dwell  time 
is  plot  in  Fig.8.  We  can  barely  see  any  feature  from  this  figure. 
After,  the  simple  gating  technique  is  used  to  filter  out  the  body 
component  from  the  original  signal,  we  begin  to  see  some  weak 
peaks  in  the  autocorrelation  function  plotted  in  Fig.9,  in  which 
it  is  still  hard  to  determine  the  period  time  tp.  The  problem  of 
the  gating  technique  is  that  it  not  only  removed  the  component 
due  to  the  target  body,  but  also  removed  some  useful  signal 
component  due  to  rotating  blades.  The  adaptive  joint  time- 
frequency  processing,  on  the  other  hand,  keeps  the  signal 
component  due  to  rotating  blades  while  removes  the  component 
due  to  the  target  body.  The  result  is  plotted  in  Fig.  10.  As  it  can 
be  seen,  multiple  sharp  peaks  occurs  at  a  fixed  time  interval  tp. 
This  interval  is  measured  to  be  0.056  seconds.  Therefore,  the 
extracted  information  indicate  that  the  rotation  rate  times  the 
number  of  the  blades  for  this  target  is  equal  to  17.7  rps. 


Fig.9.  The  autocorrelation  function  versus  dwell  time  of 
the  signal  with  body  component  gated  out. 


Fig.8.  The  autocorrelation  function  versus  dwell  time  of 
the  original  signal. 


Fig.  10.  The  autocorrelation  function  versus  dwell  time  of 
the  signal  processed  by  the  AJTF  technique. 


5.  SUMMARY 

In  this  paper  we  have  applied  the  adaptive  joint  time-frequency  concept  to  separate  two  different  scattering  mechanisms,  one 
due  to  the  target  body  and  one  due  to  the  fast  rotating  parts  on  the  target.  This  is  accomplished  by  taking  advantage  of  their 
different  (dwell  time)-(Doppler  frequency)  behaviors.  After  adaptively  searching  for  the  linear  chirp  bases  which  best 
represent  the  time-frequency  behavior  of  the  signal  and  fully  parameterizing  the  signal  with  these  basis  functions,  we  can  sort 
these  bases  into  the  signal  component  due  to  the  fast  rotating  parts  and  the  component  due  to  target  body  motion  according 
to  their  slope  and  displacement  parameters.  The  Doppler  smearing  interference  caused  by  the  fast  rotating  parts  can  be 
removed  from  the  original  ISAR  image  by  keeping  only  those  bases  associated  with  target  body  scattering.  The  algorithm  has 
been  tested  for  numerical  simulated  radar  scattering  data  for  a  helicopter  model.  The  ISAR  image  obtained  by  removing  the 
Doppler  smearing  from  the  original  image  agrees  well  with  the  image  simulated  for  the  target  without  the  rotating  blades  on. 
Furthermore,  the  algorithm  has  also  been  applied  to  eliminate  the  body  interference  and  extract  useful  information  about  the 
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blade  rotation  rate  from  the  Doppler  signal  measured  from  a  helicopter.  Comparison  with  the  simple  gating  technique  shows 
the  effectiveness  of  this  algorithm. 
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Summary 

We  present  results  of  applying  adaptive  joint  time- 
frequency  (AJTF)  processing  to  generate  focused  ISAR 
images  from  the  TIRA  radar  data  taken  in  Germany  in 
November  1997.  Our  algorithm  utilizes  the  AJTF 
technique  to  perform  Doppler  frequency  tracking.  By 
using  a  search  and  projection  technique  in  the  joint  (dwell 
time)-(Doppler  frequency)  plane,  the  prominent  point 
scatterers  are  automatically  selected  and  tracked.  The 
higher-order  translation  and  rotation  motions  are  then 
extracted  and  compensated  for  in  the  data  to  form  a 
focused  image  of  the  target.  The  motion  compensated 
images  are  compared  against  the  reference  images 
generated  by  using  the  motion  information  available  in 
the  instrumented  ARDS  data.  Furthermore,  comparison  is 
also  made  with  the  simulated  ISAR  images  of  the  air 
target  from  the  radar  signature  prediction  code  Xpatch. 
These  results  provide  an  assessment  of  the  current 
capabilities  in  using  high-resolution  ISAR  imagery  for 
non-cooperative  target  recognition. 


/.  Introduction 

Inverse  synthetic  aperture  radar  (ISAR)  imaging  is  a 
promising  method  currently  being  investigated  for  non- 
cooperative  target  identification  (NCTI).  The  basic 
framework  of  ISAR-based  NCTI  involves  the  collection 
of  multiple  frequency,  multiple  look  data  from  an 
unknown  moving  target.  The  collected  data  are  processed 
using  motion  compensation  algorithms  to  form  a  focused 
range-Doppler  image  of  the  target.  The  formed  image  is 
then  compared  against  an  existing  database  consisting  of 
either  pre-collected  measurement  or  pre-computed 
synthetic  imageries.  Some  important  questions  that  must 
be  considered  in  ISAR-based  NCTI  include:  (i)  how  well 
can  motion  compensation  processing  be  carried  out  to 
form  a  focused  image  of  the  unknown  target,  and  (ii)  what 
is  the  quality  of  the  synthetic  image  predicted  using 
electromagnetic  calculations  in  comparison  to  the 


measured  data.  In  this  paper,  we  attempt  to  address  these 
two  issues  using  the  TIRA  radar  data  taken  in  Germany. 
Three  sets  of  images  are  generated.  First,  we  form  the 
motion  compensated  imagery  using  an  algorithm  we  have 
developed  previously  [1].  It  is  based  on  a  search  and 
projection  technique  in  the  joint  (dwell  time)-(Doppler 
frequency)  plane  to  select  and  track  the  prominent  point 
scatterers.  The  higher-order  translation  and  rotation 
motions  are  then  extracted  and  compensated  for  in  the 
data  to  form  a  focused  image  of  the  target.  Second,  we 
form  the  truth  image  by  using  the  ARDS  (Advanced 
Range  Data  System)  motion  data  obtained  from  GPS 
sensors  installed  on-board  the  air  targets.  We  derive  the 
range  and  look  angle  of  the  target  relative  to  the  ground 
radar  based  on  the  ARDS  drta  ~et.  We  then  utilize  this 
information  with  interpolation  and  FFT-processing  to 
obtain  a  reference  image  of  the  target.  Third,  we  generate 
the  synthetic  image  by  using  a  CAD  model  of  the  target 
and  the  radar  prediction  code  Xpatch  [2].  Xpatch  is  based 
on  the  shooting  and  bouncing  ray  technique  [3]  and 
accounts  for  the  multiply  reflected  returns  from  a  target. 
The  three  methods  of  generating  the  images  are  detailed 
in  Sections  2,  3  and  4.  We  then  present  some  quantitative 
comparisons  between  these  images  and  draw  some 
conclusions  in  Sections  5  and  6. 


2.  Motion  Compensated  Image 

To  form  a  focused  image  from  the  raw  TIRA  data,  we 
first  carry  out  a  coarse  alignment  of  the  data  in  the  range 
dimension,  followed  by  fine  motion  compensation  in  the 
cross  range  dimension.  128  pulses  are  used  to  form  the 
image  at  each  look  angel.  During  the  coarse  alignment, 
we  correlate  the  range  profiles  among  pulse  numbers. 
The  resulting  range  shifts  between  pulses  are  then  fitted  to 
a  low-order  polynomial  to  derive  a  more  accurate  target 
range  function  versus  pulse  number.  This  function  is  used 
to  align  the  range  profiles.  Next,  fine  motion 
compensation  is  carried  out  using  an  adaptive  joint  time- 
frequency  algorithm  [l].  First,  we  select  the  range  ceil 


containing  the  brightest  scatterer.  Due  to  translation 
motion,  the  Doppler  frequency  versus  dwell  time  behavior 
of  the  point  scatterers  within  this  range  cell  is  not  constant 
in  the  joint  time-frequency  plane  (see  Fig.  1).  To  find  the 
motion  parameters,  we  project  the  response  within  the 
•range  cell  to  basis  functions  of  the  form 

/i(7)  =  expf-  +  j/'F  +y/V  +  ...)] 


Range  profiles  versus  pulse  no. 


Doppler  Frequency 


Fig.  1.  Fine  motion  compensation  is  carried  out  by  extracting 
the  Doppler  frequency  versus  dwell  time  behavior  of 
the  strong  point  scatterer  in  the  signal. 

where  t  denotes  dwell  time  (or  pulse  number).  By 
varying  ff  and  /  and  searching  for  the  maximum 
projection  value,  the  basis  function  which  best  resembles 
the  joint  (dwell  time)-(Doppler  frequency)  behavior  of  the 
strongest  point  scatterer  in  the  range  cell  can  be  found. 
For  TIRA  data,  we  use  polynomials  up  to  the  third  order 
in  the  phase  function  of  the  basis.  The  search  of  the  first 
order  coefficient  /  can  be  accomplished  by  using  an  FFT 
algorithm,  while  a  brute-force  search  is  needed  to  find  the 
higher  order  coefficients /  and  /  .  Since  the  higher  order 
phase  terms  in  the  basis  correspond  to  phase  errors  cause 
by  the  motion,  we  can  multiply  the  conjugate  of  this  basis 
function  to  the  original  signal  to  eliminate  the  translation 
motion  from  it.  This  algorithm  can  also  be  extended  to 


multiple  range  ceils  to  correct  for  higher  order  rotation 
motions,  although  we  did  not  find  them  to  be  significant. 
Shown  in  Figs.  3(a),  4(a)  and  5(a)  are  the  motion 
compensated  images  from  the  TIRA  data  at  three  different 
look  angles  processed  using  the  above  motion 
compensation  algorithm. 


3.  Reference  Image  Based  on  ARDS  Motion  Data 

Several  targets  within  the  TIRA  data  set  carried  on¬ 
board  GPS  sensors  to  record  the  position  and  motion  of 
the  aircraft.  The  data  were  processed  and  written  in  the 
ARDS  file.  We  form  images  based  on  the  ARDS  data  set 
and  use  them  as  a  reference  to  evaluate  the  quality  of  the 
motion  compensated  images  described  in  the  last  section. 
The  ARDS  data  set  provides  sensor-derived  range  and 
orientation  information  (yaw,  pitch  and  roll)  of  the  target. 
We  use  the  range  information  to  carry  out  range 
alignment  after  first  synchronizing  the  GPS  time  and  the 
radar  collection  time.  The  range  data  are  then  fitted  to  a 
low-order  polynomial  in  the  same  manner  as  the  coarse 
motion  compensation  procedure.  The  resulting  range 


Fig.  2.  Elevation  and  azimuth  look  angels  on  a  TIRA 
target  during  its  flight. 

function  is  used  to  align  the  range  profiles  over  pulse 
numbers.  To  carry  out  Doppler  alignment,  the  target 
orientation  (yaw,  pitch  and  roll)  is  first  transformed  into 
the  radar  look  angle  on  target  in  terms  of  elevation  and 
azimuth.  Fig.  '2  shows  the  elevation  and  azimuth  look 
angles  on  a  target  during  the  dwell  history  of  11.7 
minutes.  It  can  be  seen  that  in  this  case  the  elevation 
angle  is  near  zero  during  the  entire  flight.  Given  the 
relative  constant  elevation  angle  over  the  entire  data 
collection  time,  we  perform  Doppler  focusing  by 
assuming  only  azimuth  variation.  Since  the  imaging 
plane  is  the  zero-elevation  plane  under  this  condition,  the 
expected  2D  ISAR  images  will  be  the  top  view  of  the 
aircraft.  The  azimuth  versus  dwell  time  information  is 


first  fitted  to  a  low-order  polynomial  and  the  resulting 
function  is  used  to  reformat  the  data  from  uniform 
sampling  in  dwell  time  to  uniform  sampling  in  azimuth 
angle.  FFT  processing  is  then  applied  to  form  the  image. 
Figs.  3(b),  4(b)  and  5(b)  show  the  reference  images 
-derived  from  the  ARDS  motion  data  for  the  same  target  at 
the  same  looks  as  those  shown  in  Figs.  3(a),  4(a)  and  5(a). 
128  pulses  are  used  to  form  the  images. 


4.  Synthetic  Image  from  Xpatch  Simulation 

We  simulate  the  synthetic  imagery  based  on  a  CAD 
model  of  the  target  and  using  the  radar  simulation  code 
Xpatch  [2].  The  CAD  model  consists  of  a  facetized 
version  of  the  target.  In  the  Xpatch  simulation, 
geometrical  optics  rays  are  shot  from  the  incident  look 
angle  and  all  the  multiple  reflections  are  tracked  until  the 
rays  exit  the  target.  The  image  is  formed  by  updating  the 
ISAR  image  plane  one  ray  at  a  time  using  its  ray-spread 
function.  A  fast  FFT-based  algorithm  is  next  used  to 
accelerate  the  ray  update  time.  This  fast  image  formation 
algorithm  is  based  on  a  single  ray  trace  and  is  orders  of 
magnitude  faster  than  the  conventional  multi-frequency, 
multi-aspect  ray  calculation.  It  has  been  found  from 
validation  studies  that  this  method  leads  to  a  more 
focused  image  than  the  conventional  frequency-aspect 
approach.  However,  the  key  image  features  are  fairly 
well  predicted  by  this  algorithm.  The  detailed  description 
of  the  algorithm  can  be  found  in  [4]  and  [5].  The  typical 
simulation  time  per  image  (carried  out  at  a  ray  density  of 
10  rays/wavelength)  for  a  faceted  aircraft  model  at  Ku 
band  is  approximately  30  minutes  on  an  SGI  02 
workstation.  The  simulation  time  is  dominated  by  the  ray 
tracing  operation.  The  single  diffraction  contribution 
from  edges  can  also  be  computed  in  the  same  manner 
using  Xpatch.  However,  we  did  not  include  it  in  our 
calculations.  Figs.  3(c),  4(c)  and  5(c)  show  the  results  of 
the  Xpatch  simulation  for  the  same  target  at  the  same  set 
look  angles  as  the  measured  images.  Both  the  down 
range  and  cross  range  resolutions  in  the  simulation  are 
adjusted  to  correspond  to  the  measured  data. 


5.  Image  Comparisons 

We  have  processed  the  images  of  a  TIRA  target  over 
its  flight  duration  of  11.7  minutes.  Three  sets  of 
representative  examples  are  shown  here.  Figs.  3(a),  (b) 
and  (c)  are  the  images  of  the  target  at  azimuth=146°  (0° 
being  nose-on)  generated  using,  respectively,  motion 
compensation,  ARDS  sensor  information  and  Xpatch 
simulation.  The  dynamic  range  of  the  displayed  images  is 
55  dB.  Figs.  4(a),  (b)  and  (c)  are  the  images  of  the  target 
at  azimuth=-56°  obtained  by  the  three  methods.  Figs.  5(a), 
(b)  and  (c)  are  the  images  of  the  target  at  azimuth=30°. 
The  look  angle  information  was  deciphered  from  the 
ARDS  data.  We  make  the  following  observations: 


(1)  The  motion  compensated  images  and  the  ARDS- 
derived  reference  images  appear  to  be  in  good  agreement 
over  most  angles.  Surprisingly,  the  motion  compensated 
images  are  better  focused  than  the  reference  images.  We 
attribute  this  to  the  finite  accuracy  of  the  GPS  sensors, 
which  is  estimated  to  be  about  6  cm.  Since  this 
inaccuracy  amounts  to  several  wavelengths  at  Ku  band, 
we  believe  the  motion  data  did  not  provide  sufficient 
accuracy  for  forming  a  well-focused  image.  On  the  other 
hand,  we  found  that  the  performance  of  our  motion 
compensation  algorithm  depends  on  the  availability  of  a 
clear  point  scatterer  in  the  selected  range  cell.  A  general 
criterion  should  be  further  developed  to  automatically 
select  the  range  cell  that  contains  a  strong,  well-isolated 
point  scatterer.  Shown  in  Fig.  6(a)  is  the  correlation 
coefficient  between  the  two  sets  of  images  over  40  looks 
that  almost  cover  the  entire  azimuth  scan  on  the  target. 
The  correlation  coefficient  is  normalized  between  zero 
and  one  and  is  computed  based  on  the  image  intensity 
with  a  power  transform  of  0.5.  We  observe  that  the 
correlation  coefficient  is  on  the  order  of  0.9  over  the 
azimuth  sweep. 

(2)  The  Xpatch  simulation  and  the  ARDS-derived 
reference  images  show  fairly  good  agreement  in  the 
prominent  target  features.  As  expected,  the  Xpatch 
images  are  more  focused  and  do  not  exhibit  the  diffused 
characteristics  of  the  measurement  data.  This  observation 
is  consistent  with  other  data  sets  we  have  studied  in  the 
past.  We  also  observe  that  in  the  frontal  look  region  of 
the  aircraft,  the  measured  images  contain  strong  jet  engine 
modulation  (JEM)  lines  as  shown  in  Fig.  5.  This  effect  is 
not  predicted  in  the  Xpatch  simulation.  Fig.  6(b)  shows 
the  correlation  coefficient  between  the  two  sets  of  images 
as  a  function  of  azimuth  angle.  The  correlation  between 
the  synthetic  data  and  the  measured  reference  data  is  on 
the  order  of  0.4  to  0.5.  Near  the  frontal  region  (±50°),  the 
correlation  is  lower  due  to  JEM  effects  in  the  measured 
data.  We  have  recently  developed  a  methodology  to 
predict  JEM  lines  using  Xpatch  [6].  However,  this  effect 
is  strongly  dependent  on  the  engine  spin  rate  relative  to 
the  radar  pulse  repetition  frequency  and  an  accurate 
prediction  is  difficult. 

6.  Conclusions 

We  have  reported  some  preliminary  results  from 
processing  the  TIRA  measurement  data.  We  have  applied 
the  adaptive  joint  time-frequency  motion  compensation. to 
form  ISAR  imageries  of  TIRA  targets.  The  results 
compared  favorably  with  the  reference  images  formed 
using  the  ARDS  motion  data.  It  was  found  that  the 
quality  of  the  reference  image  is  slightly  inferior  to  the 
motion  compensated  image  due  to  the  accuracy  limit  of 
the  GPS  sensors.  While  the  motion  compensation 
algorithm  performed  quite  adequately  in  this  instance,  the 
target  under  consideration  did  not  exhibit  significant 
higher-order  rotation  motion.  Faster  maneuvering  targets 
may  pose  more  of  a  challenge  to  the  motion  compensation 


algorithm.  Also,  we  have  arbitrarily  chosen  128  pulses  to 
form  the  images.  We  will  investigate  the  use  of  longer 
dwell  time  in  the  image  formation  process  to  achieve 
better  cross  range  resolution  and  extract  absolute  cross 
range  scaling. 

We  have  also  generated  the  synthetic  images  of  the 
target  and  compared  the  synthetic  results  to  the  reference 
images  from  the  TIRA  data.  It  was  found  that  a  fairly 
good  agreement  was  achieved  using  Xpatch  in  predicting 
the  prominent  target  features.  A  higher-fidelity  CAD 
model  should  further  improve  the  quality  of  the  synthetic 
data.  Within  the  frontal  sector  of  the  target,  strong  JEM 
lines  in  the  measured  images  significantly  corrupted  the 
geometrical  features  on  the  target.  We  are  pursuing 
algorithms  for  removing  such  artifacts  from  the  measured 
data  to  unmask  the  geometrical  features  of  the  target  [7]. 
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Fig.  3.  A  TIRA  target  at  azimuth=146°. 

(a)  Motion  compensated  image  using  AJTF  processing. 

(b)  Reference  image  based  on  ARDS  motion  data. 

(c)  Synthetic  image  from  Xpatch  simulation. 
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Fig.  4.  A  TIRA  target  at  azimuth=-56°. 

(a)  Motion  compensated  image  using  AJTF  processing. 

(b)  Reference  image  based  on  ARDS  motion  data. 

(c)  Synthetic  image  from  Xpatch  simulation. 
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Fig.  6.  Correlation  coefficients  versus  azimuth  look  angle  between: 

(a)  Motion  compensated  images  and  ARDS-derived  reference  images. 

(b)  Xpatch  simulated  images  and  ARDS-derived  reference  images. 
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l.Introduction 

Synthetic  Aperture  Radar  (SAR)  images  of  ground  targets  generally  consist 
of  target  features  and  clutters  from  background  scattering  [1].  In  automatic 
target  recognition  (ATR)  applications,  it  is  desirable  to  remove  the  clutter  from 
the  actual  target  images  before  ATR  processing.  The  standard  way  to  suppress 
clutter  is  to  apply  an  appropriate  threshold  level  to  the  whole  SAR  image  to 
remove  the  clutter.  However,  this  approach  assumes  that  the  target  signal-to- 
clutter  ratio  (SCR)  is  large  enough,  otherwise  results  in  either  some  target 
feature  loss  or  remnant  clutter  residue.  In  this  work,  we  set  out  to  develop  a 
decluttering  algorithm  to  automatically  extract  the  target  image  from  a  SAR 
image  by  maximizing  SCR  using  the  adaptive  wavelet  packet  transform 
(AWPT)  [2].  The  wavelet  packet  basis  is  the  generalization  of  the  conventional 
wavelet  basis  [3]  and  has  been  applied  for  image  compression  [4]  and  moment 
matrix  sparsification  [2]. 

Our  approach  is  to  transform  the  SAR  image  to  a  new  domain  using  an 
appropriately  chosen  set  of  wavelet  packet  basis.  Since  a  typical  target  image 
usually  consists  of  point  scatterers  and  region  fixtures,  which  are  respectively 
high  frequency  and  low  frequency  signals,  the  multi-scale  wavelet  basis  is  well 
suited  to  focus  the  target  image.  Clutter  image,  on  the  other  hand,  is  statistically 
uncorrelated  from  pixel  to  pixel.  Therefore,  the  transformed  clutter  image  under 
the  same  set  of  bases  should  remain  uncorrelated.  We  therefore  expect  that  the 
SCR  can  be  increased  in  an  appropriately  chosen  set  of  basis.  The  cost  function 
of  our  AWPT  algorithm  is  chosen  to  describe  how  well  the  target  signal  is 
focused  in  the  transform  domain.  An  efficient  basis  search  algorithm  is 
implemented  to  find  the  best  wavelet  packet  basis.  Our  algorithm  is  tested  using 
the  MSTAR  SAR  data  set  [5]  and  show  that  an  improved  SCR  can  be  achieved 
using  the  algorithm. 

2.  SAR  Image  Representation  with  Wavelet  Packet  Basis 

A  discrete  SAR  image  s(m,  n)  can  be  represented  as: 

s(m,n)  =  t(m,n)  +  c(m,n )  0<n,m<  N  (1) 

where  t(m,  n)  and  c(m,  n)  denote  the  target  image  and  the  clutter  in  the  SAR 
image,  respectively.  We  define  a  set  of  orthogonal  and  complete  2-D  wavelet 
packet  basis  functions: 

{U'pq(kJ)  0<  j  <J,  0<  p,q  <2' ,  0<k,I  <  N2~') 


(2) 


where  j  denotes  the  scale  index,  and  J  =  Log2(N);  p,  q  are  the  modulation 
indices;  and  k,  1  are  the  position  indices.  The  2-D  wavelet  packet  basis  function 
can  be  generated  from  the  product  between  two  1-D  wavelet  packet  bases: 

UJp,g  (k,l)= WJP  (k)yfJv  (/)  (3) 

\|/  is  a  1-D  wavelet  packet  basis  that  can  be  generated  from  the  scaling  function 
and  the  basic  wavelet  function  using  the  “2-scale  equation”  [2-4]. 

The  transformation  of  a  SAR  image  s(m,  n)  using  the  wavelet  packet  basis  is: 
SJM(k,l)  =  ZXs(m,n)UJpJk-  ym,l-2Jn) 

m  n 

=  ZZt(m,n)Up  c/(k- ?  m,l  -  21  n)  +  Z Z t(m,n)lfp  q(k  -  21  m,l  -  21  n) 

m  n  mn 


=  fUk,i)  +  cUk,i) 


(4) 

where  T  and  C  are  the  transform  coefficients  of  the  target  image  and  the  clutters 
in  the  image  with  the  wavelet  packet  basis. 

If  we  define  SCR  of  a  SAR  image  as  the  ratio  of  average  target  amplitude  to 
the  standard  deviation  of  the  clutter,  the  SCR  of  the  image  before  and  after  the 

transform  are  1 1  <j(c)andT  /<r(C)  respectively.  The  clutters  function  c(m,  n) 
denotes  the  reflectivity  of  different  resolution  cells  in  the  highlighted 
background  areas,  and  are  independently  and  identically  distributed  [6].  It  can 

be  shown  that  their  wavelet  packet  transform  coefficients  C  are  still 
uncorrelated  [7],  and  the  standard  deviation  of  clutter  does  not  change  after  the 
basis  transform.  But  the  signals  from  target  area  are  related  to  themselves,  it  is 
possible  to  increase  target  amplitude  with  a  basis  transformation.  Therefore, 
with  a  specific  SAR  image  we  need  to  find  the  best  wavelet  packet  basis  to 
maximize  the  transformed  target  image  magnitude  adaptively. 


3.  Adaptive  Search  Procedure  and  Implementation 

To  find  the  best  wavelet  packet  basis  and  implement  the  basis  transform,  we 
need  to  define  a  cost  function  to  describe  how  well  the  transformed  target  signal 
is  concentrated  with  a  wavelet  packet  basis.  The  best  basis  is  the  one  that  makes 
the  transformed  signal  have  the  minimum  cost.  The  most  commonly  used  cost 
function  is  the  entropy  function.  Because  of  the  complexity  of  evaluating 
entropy  function,  we  use  energy  concentration  function  as  cost  function  in  this 
application.  For  a  transformed  SAR  image  it  is  defined  as: 

ip 


Cost  =  zz 

j.p,q  k,l 


si. 


0<  P<  2 


(5) 


Because  there  are  many  possible  wavelet  packet  bases  in  (2)  for  the  transform,  it 
is  impractical  to  try  each  of  them  to  find  the  best  basis.  An  effective  Quad-tree 
decomposition  algorithm  was  generalized  from  that  proposed  in  [8]  to  find  the 
best  2-D  wavelet  packet  basis  with  regard  to  an  additive  cost  function.  The 
algorithm  decomposes  the  original  image  using  2-channel  filtering  with  a  pair  of 
Quadrature  filters  stage  by  stage  from  spatial  domain  to  spectrum  domain.  With 
cost  labeled  at  every  stage,  the  decomposition  tree  is  searched  from  the  last  stage 
to  the  first  stage.  The  best  decomposition  tree  can  be  found  after  the  search 


process  with  total  computational  complexity  of  about  0(N2logN)  for  the 
algorithm.  A  scale-dependent  threshold  is  applied  to  the  transformed  image. 
Because  there  is  a  weak  correlation  between  clutter  samples,  we  increase 
threshold  level  slightly  as  the  scale  increases.  The  threshold  level  is  chosen  as: 

Thj  =  Kacy[JTJ  (6) 

where  j  is  the  scale  index,  J  is  the  maximum  scale,  cc  is  the  clutter  standard 
deviation,  and  K  is  a  constant.  With  the  thresholding  processing,  most  clutters 
are  removed  in  the  wavelet  packet  basis  domain,  and  the  image  is  inverse- 
transformed  back  to  spatial  domain  to  restore  the  original  target  image  using  the 
same  tree  for  decomposition.  Although  SCR  is  much  improved  in  the  restored 
image,  it  is  impossible  to  remove  the  clutter  completely  through  that  processing. 
We  apply  a  very  small  second  threshold  to  the  restored  image,  and  then 
clustering  processing  to  get  rid  of  all  clutter  residues. 

4.  Test  Results 

To  test  the  effectiveness  of  the  processing  scheme,  we  apply  it  to  the  MSTAR 
SAR  image  data.  Fig.  1  shows  one  of  MSTAR  images  in  which  target  is  a 
ground  vehicle  with  vegetation  clutters.  There  are  several  strong  point  scatters  in 
the  front  of  the  vehicle,  but  the  reflection  from  the  back  part  is  relatively  weak. 
With  direct  thresholding  processing,  as  demonstrated  in  Fig.2,  some  crucial 
parts  of  the  target  are  lost.  However,  if  we  apply  the  AWPT  algorithm  to  the 
same  image,  the  restored  image  is  shown  in  Fig.  3.  All  key  features  of  the 
original  image  are  kept  through  the  processing.  We  chose  Daubechies  filter 
with  order  of  6  as  wavelet  filter,  the  first  threshold  parameter  is  1  and  the  second 
threshold  is  0.05ac.  In  both  processing  methods,  there  is  some  target  information 
loss.  Fig.  4  shows  the  signal-to-clutter  ration  with  average  target  amplitude  loss 
for  the  two  processing  methods.  It  is  observed  that  for  a  fixed  target  image  loss 
the  new  AWPT  method  always  achieves  a  higher  SCR  value  than  the  traditional 
direct  thresholding  processing.  Similar  results  are  obtained  when  the  algorithm 
is  applied  to  other  MSTAR  data. 
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Fig.3.  Clutter  rejection  using  AWPT  Fig.  4.  SCR  vs.  target  image  loss  for 
method  the  two  processing  methods. 


